{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Approximating Derivatives by the Method of Undetermined Coefficients" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Last revised Monday March 1,** adding several more worked examples." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**References:**\n", "\n", "- Section 5.1 *Numerical Differentiation* of [Sauer](../references.html#Sauer)\n", "- Section 4.1 *Numerical Differentiation* of [Burden&Faires](../references.html#Burden-Faires)\n", "- Section 4.2 *Estimating Derivatives and Richardson Extrapolation* of [Chenney&Kincaid](../references.html#Chenney-Kincaid)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have seen several formulas for approximating a derivative $Df(x)$ or higher derivative\n", "$D^k f(x)$ in terms of several values of the function $f$, such as" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$ Df(x) \\approx D_hf(x) := \\frac{f(x+h) - f(x)}{h} $$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$ D^2f(x) \\approx \\delta^2 f(x) := \\frac{f(x-h) - 2f(x) + f(x+h)}{h^2}. $$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the first case we get an error formula" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$D_hf(x) = \\frac{f(x+h) - f(x)}{h} = Df(x) + \\frac{D^2f(\\xi)}{2} h, \\quad \\xi \\text{ between $x$ and $x+h$}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "by inserting the Taylor formula $f(x+h) = f(x) + Df(x) + \\frac12 {D^2f(\\xi)}$,\n", "and thus an \"order of accuracy formula\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$D_hf(x) = \\frac{f(x+h) - f(x)}{h} = Df(x) + O(h)$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "so that the error is" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$ D_hf(x) - Df(x) = \\frac12 {D^2f(\\xi)} = O(h). $$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These are linear combinations of values of $f$ at various points, with the denominator scaling with the k-th power of the mode spacing scale $h$, which makes sense given the linearity of derivatives and the way that the k-th derivative scales when one rescales $f(x)$ to $f(ck)$.\n", "\n", "Thus we will make the *Ansatz* that the k-th derivative $D^k f(x)$ can be approximated using values at th $r-l+1$ equally spaced points" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$x + lh, x + (l+1)h, \\dots x + rh$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where the integers $l$ and $r$ can be negative, positive or zero:\n", "the assumed form then is" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$D^k f(x) \\approx D_h^k f(x) = \\frac{C_l f(x + lh) + C_{l+1} f(x + (l+1)h) + \\cdots + C_r f(x + rh)}{h^k} + O(h^p)$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(The reason for the power $k$ in the denominator will be seen soon.)\n", "\n", "So we seek to determine the values of the initially undetermined coefficients $C_i$,\n", "by the criterion of giving an error $O(h^p)$ with the highest possible order $p$.\n", "With $r-l+1$ coefficients to choose, we generally get $p = r - l + 1 - k$, but with symmetry $l = -r$ and $k$ even we get one better, $p = r - l + 2 - k$, because the order $p$ must then be even.\n", "Thus we need the number of points $r-l+1$ to be more than $k$: for example, at least two for a first derivative as already seen." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Example 1**\n", "\n", "The basic forward difference approximation\n", "\n", "$$ Df(x) = \\frac{f(x + h) - f(x)}{h} + O(h) $$\n", "\n", "has $k=1$, $l=0$, $r=1$, $p=1$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Example 2**\n", "\n", "A three-point one-sided difference approximation of the first derivative ($k=1$) can be sought with $l=0$, $r=2$, as\n", "\n", "$$Df(x) = \\frac{C_{0} f(x) + C_1 f(x + h) + C_2 f(x + 2h)}{h} + O(h^p)$$\n", "\n", "and the most accurate choice is $C_0 = -3/2$, $C_1 = 2$, $C_2 = -1/2$, again of second order, which is exactly\n", "$p = r - l + 1 - k$, with no \"symmetry bonus\":\n", "\n", "$$Df(x) \\approx \\frac{-3 f(x) + 4 f(x + h) - f(x + 2h)}{2 h} + O(h^2).$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Example 2, continued**\n", "\n", "One can use Taylor's Theorem to *check* an approximation like this, and also get information about its accuracy.\n", "To do this, insert a Taylor series formula with center $x$, like\n", "\n", "$$f(x+h) = f(x) + Df(x) h + \\frac{D^2f(x)}{2} h^2 + \\frac{D^3f(x)}{6} h^3 + \\cdots$$\n", "\n", "If you are not sure how accurst the result is, you might need to initially be vague about how may terms are needed,\n", "so I will do it that way and then go back and be more specific once we know more.\n", "\n", "A series for $f(x+2h)$ is also needed:\n", "\n", "$$\\begin{split}\n", "f(x+2h) &= f(x) + Df(x) (2 h) + \\frac{D^2f(x)}{2} (2 h)^2 + \\frac{D^3f(x)}{6} (2 h)^3 + \\cdots\n", "\\\\&= f(x) + 2 Df(x) h + \\frac{D^2f(x)}{2} 4 h^2 + \\frac{D^3f(x)}{6} 8 h^3 + \\cdots\n", "\\\\&= f(x) + 2 Df(x) h + 2 D^2f(x) h^2 + \\frac{4 D^3f(x)}{3} h^3 + \\cdots\n", "\\end{split}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Insert these into the above three-point formula, and see how close it is to the exact derivative:\n", "\n", "$$\\begin{split}\n", "&\\frac{-3 f(x) + 4 f(x + h) - f(x + 2h)}{2 h}\n", "\\\\&= \\frac{-3 f(x) + 4 [f(x) + Df(x) h + \\frac{D^2f(x)}{2} h^2 + \\frac{D^3f(x)}{6} h^3 + \\cdots]\n", "- [f(x) + 2 Df(x) h + 2 D^2f(x) h^2 + \\frac{4 D^3f(x)}{3} h^3 + \\cdots]}{2 h}\n", "\\end{split}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now gather terms with the same power of $h$ (which is also gathering terms with the same order of derivative):\n", "\n", "$$\\begin{split}\n", "\\frac{-3 f(x) + 4 f(x + h) - f(x + 2h)}{2 h}\n", "&= f(x)\\frac{-3 + 4 -1}{2h} + Df(x)\\frac{4 -2}{2} + D^2f(x)\\frac{4/4 - 2/2}h + D^3f(x)\\frac{4/12 - 4/6}h^2 + \\cdots\n", "\\\\&= Df(x) - \\frac{D^3 f(x)}{3} h^2 + \\cdots\n", "\\end{split}$$\n", "\n", "and it is clear that the omitted terms have higher power of $h$: $h^3$ and up.\n", "That is, they are $O(h^3)$, or more conveniently $o(h^2)$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Thus we have confirmed that the error in this approximation is\n", "\n", "$$Df(x) - \\frac{-3 f(x) + 4 f(x + h) - f(x + 2h)}{2 h} = \\frac{D^3 f(x)}{3} h^2 + o(h^2) = O(h^2).$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Example 3**\n", "\n", "A three-point centered difference approximation of $D^2 f(x)$ [$k=2$] has $l = -1$, $r = 1$ and so\n", "\n", "$$D^2f(x) \\approx \\frac{C_{-1} f(x - h) + C_{0} f(x) + C_1 f(x + h)}{h^2}$$\n", "\n", "and it can be found (as discussed below) that the coefficients $C_{-1} = C_1 = 1$, $C_0 = -2$ give the highest order error: $p=2$; one better than $p = r - l + 1 - k = 1$ due to symmetry:\n", "\n", "$$D^2f(x) = \\frac{f(x - h) -2 f(x) + f(x + h)}{h^2} + O(h^2).$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Method 1: use Taylor polynomials in $h$ of degree $p+k-1$, so with error terms $O(h^{p+k})$\n", "\n", "Each of the terms $f(x + ih)$ in the above formula for the approximation $D_h^k f(x)$ of the $k$-th derivative $D_k f(x)$\n", "can be expanded with the Taylor Formula up to order $p+k$," ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$f(x + ih) = f(x) + (ih)Df(x) + (ih)^2/2 D^2f(x) + \\cdots + (ih)^j/j! D^jf(x) + \\cdots + (ih)^{p+k}/(p+k)! D^{p+k}f(x) + o(h^{p+k})$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then these can be rearranged, putting the terms with the same derivative $D^j f(x)$ together —\n", "all of which have the same factor $h^j$ in the numeriator, and so the same factor $h^{j-p}$ overall:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\\begin{split}\n", "D_h^k f(x)\n", "&= (C_l + \\cdots + C_r)f(x)h^{-k}\\\\\n", "&+ (l C_l + \\cdots + r C_r)Df(x)h^{1-k}\\\\\n", "&+ (l^2 C_l + \\cdots + r^2 C_r)D^2f(x)h^{2-k}\\\\\n", "& \\vdots\\\\\n", "&+ (l^j C_l + \\cdots + r^j C_r)D^jf(x)h^{j-k}\\\\\n", "& \\vdots\\\\\n", "&+ (l^{p+k} C_l + \\cdots + {p+k}^j C_r)D^{p+k}f(x)h^{p}\\\\\n", "&+ o(h_p)\n", "\\end{split}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The final \"small\" term $o(h^{p})$ comes from the terms $o(h^{p+k})$ in each Taylor's formula term, each divided by $h^k$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We want this whole thing to be approximately $D^k f(x)$,\n", "and the strategy is to match the coefficients of the derivatives:\n", "\n", "1. Matching the coefficients of $D_h^k f(x)$," ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$(l^k C_l + \\cdots + r^k C_r) D^k f(x)h^{k-k} = (l^k C_l + \\cdots + r^k C_r) D^k f(x) = D^k f(x)$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "so" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$l^k C_l + \\cdots + r^k C_r = 1$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "2. On the other hand, there should be no term with factor $f(x)h^{-k}$, so" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$C_l + \\cdots + C_r = 0$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "More generally, for any $j$ other than $k$ the coefficients should vanish, so" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$l^j C_l + \\cdots + r^j C_r = 0, \\quad 0 \\leq j \\leq p+k \\text{ except for } j = k$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This last line gives $p+k$ linear equations in the $p+k+1$ coefficients $C_1 , \\dots, C_{p+k}$,\n", "and then the previous equation gives us a total of $p+k+1$ equations — as needed for the existence of a unique solution." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{split}\n", "C_l + \\cdots + C_r &= 0\\\\\n", "l^j C_l + \\cdots + r^jC_r &= 0, j \\neq k\\\\\n", "l^k C_l + \\cdots + r^kC_r &= 1\\\\\n", "\\end{split}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And indeed it can be verified that the resulting matrix for this system of equations is non-singular, and so there is a unique solution for the coefficients $C_l \\dots C_r$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### A) Derive the formula in Example 2.\n", "\n", "Do this by setting up the three equations as above for the coefficients $C_0$, $C_1$ and $C_2$, and solving them.\n", "Do this \"by hand\", to get exact fractions as the answers;\n", "us the two Taylor serei formulas, but now tak advanta of what we saw above, whi cis that the error stsrts at the terms in\n", "$D^3f(x)$, so use the forms\n", "\n", "$$f(x+h) = f(x) + Df(x) h + \\frac{D^2f(x)}{2} h^2 + \\frac{D^3f(x)}{6} h^3 + O(h^4)$$\n", "\n", "and\n", "\n", "$$f(x+2h) = f(x) + 2 Df(x) h + 2 D^2f(x) h^2 + \\frac{4 D^3f(x)}{3} h^3 + O(h^4)$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### B) Verify the result in Example 3.\n", "\n", "Again, do this by hand, and exploit the symmetry.\n", "Note that it works a bit better than expected, due to the symmetry." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Degree of Precision and testing with monomials\n", "\n", "This concept relates to a simpler way of determining the coefficients.\n", "\n", "The **degree of precision** of an approximation formula (of a derivative or integral) is the highest degree $d$ such that the formula is exact for all polynomials of degree up to $d$.\n", "For example it can be checked that in the examples above, the degrees of precision are 1, 2, and 3 respectively.\n", "All three conform to a general pattern:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Theorem:** $d = p + k - 1$, so the typical case with no \"symmetry bonus\" $d = r - l$\n", "\n", "This is confirmed by the above derivation: for $f$ any polynomial of degree $p+k-1$ or less,\n", "the Taylor polynomials of degree at most $p+k-1$ used there have no error.\n", "\n", "Thus for example, the minimal symmetric aproximation of a fourth derivative, which must have even order $p=2$, will have degree of precision 5." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Method 2: use monomials of degree up to $p+k-1$\n", "\n", "From the above degree of precision result, one can determine the coefficients by requiring degree of precision $p+k-1$, and for this it is enough to require exactness for each of the simple monomial functions $1$, $x$, $x^2$, and so on up to $x^{p+k-1}$.\n", "\n", "Also, this only needs to be tested at $x=0$, since \"translating\" the variables does not effect the result.\n", "\n", "This is probably the simplest method in practice." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Example 4 (Example 2 revisited)**\n", "\n", "The goal is to get exactness in\n", "\n", "$$\\frac{C_{0} f(x) + C_1 f(x + h) + C_2 f(x + 2h)}{h} = Df(x)$$\n", "\n", "for the monomials $f(x) = 1$, $f(x) = x$, and so on, to the highest power possible,\n", "and this only needs to be checked at $x=0$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, $f(x) = 1$, so $Df(0) = 0$:\n", "\n", "$$\\frac{C_{0} \\times 1 + C_1 \\times 1 + C_2 \\times 1}{h} = 0,$$\n", "\n", "so\n", "\n", "$$C_{0} + C_1 + C_2 = 0$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, $f(x) = x$, so $Df(0) = 1$:\n", "\n", "$$\\frac{C_{0} f(0) + C_1 f(h) + C_2 f(2h)}{h} = \\frac{C_{0} 0 + C_1 h + C_2 2h}{h} = C_1 + 2C_2 = 1$$\n", "\n", "so\n", "\n", "$$C_1 + 2 C_2 = 1$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We need at least three equations for the three unknown coefficients, so continue with $f(x) = x^2$, $Df(0) = 0$:\n", "\n", "$$\\frac{C_{0} f(0) + C_1 f(h) + C_2 f(2h)}{h} = \\frac{C_{0} 0 + C_1 h^2 + C_2 (2h)^2}{h} = (C_1 + 4 C_2)h = 0$$\n", "\n", "so\n", "\n", "$$C_1 + 4 C_2 = 0$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can solve these by elimination; for example:\n", "- The last equation gives $C_1 = -4C_2$\n", "- The previous one then gives $-4C_2 + 2C_2 = 1$, so $C_2 = -1/2$ and thus $C_1 = -4C_2 = 2$.\n", "- The first equation then gives $C_0 = -C_1 - C_2 = -3/2$\n", "all as claimed above." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So far the degree of precision has been shown to be at least 2.\n", "In some cases it is better, so let us check by looking at $f(x) = x^3$:\n", "\n", "$Df(x) = 0$, whereas\n", "\n", "$$\n", "\\frac{-3 f(x) + 4 f(x + h) - f(x + 2h)}{2 h}\n", "= \\frac{-3 0 + 4 h^3 - (2h)^3}{2 h}\n", "= \\frac{ -2 h^3}{2 h}\n", "= -4 h^2, \\neq 0\n", "$$\n", "\n", "So, no luck this time (that typically requires some symmetry),\n", "but this calculation does indicate in a relatively simple way that the error is $O(h^2)$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you want to verify more rigorously the order of accuracy of a formula devised by this method,\n", "one can use the \"checking\" procedure with Taylor polynomials and their error terms as done in Example 2 above." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 2: like Exercise 1, but using Method 2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### A)\n", "Verify the result in Example 2, this time by Method 2.\n", "\n", "That is, impose the condition of giving the exact value for the derivative at $x=0$ for the monomial $f(x) = 1$,\n", "then the same for $f(x) = x$, and so on until there are enough equations to determine a unique solution for the coefficients." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### B)\n", "Verify the result in Example 3, by Method 2." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "This work is licensed under [Creative Commons Attribution-ShareAlike 4.0 International](https://creativecommons.org/licenses/by-sa/4.0/)" ] } ], "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.5" } }, "nbformat": 4, "nbformat_minor": 4 }