8.1 Numerical differentiation
The definition \(f'(x_0) = \lim _{h\rightarrow 0}\frac {f(x_0+h)-f(x_0)}{h}\) suggests the forward-difference formula \[f'(x_0)\approx \frac {f(x_0+h)-f(x_0)}{h}\] (the backward-difference formula when \(h<0\)). By Taylor’s Theorem, \(f(x_0+h) = f(x_0)+hf'(x_0)+\frac {h^2}2f''(\xi )\), so \[f'(x_0) = \frac {f(x_0+h)-f(x_0)}{h}-\frac h2f''(\xi ):\] the error is \(O(h)\).
Formulas from interpolation
More accurate formulas come from differentiating the Lagrange polynomial through \(n+1\) points. Differentiating the identity of Theorem 6.1.3 and evaluating at a node \(x_j\) gives \[f'(x_j) = \sum _{k=0}^nf(x_k)L'_{n,k}(x_j)+\frac {f^{(n+1)}(\xi (x_j))}{(n+1)!}\prod _{\substack {k = 0\\ k\neq j}}^n(x_j-x_k),\] because the term involving \(\frac {d}{dx}f^{(n+1)}(\xi (x))\) is multiplied by the nodal polynomial, which vanishes at \(x_j\). With three equally spaced points this gives the following formulas, with \(\xi \) in the interval spanned by the points used.
Theorem 8.1.1. Three-point endpoint formula: \[f'(x_0) = \frac {1}{2h}\bigl [-3f(x_0)+4f(x_0+h)-f(x_0+2h)\bigr ]+\frac {h^2}{3}f'''(\xi _0) .\] Three-point midpoint formula: \[f'(x_0) = \frac {1}{2h}\bigl [f(x_0+h)-f(x_0-h)\bigr ]-\frac {h^2}{6}f'''(\xi _1) .\]
Proof. For the midpoint formula, expand to third order about \(x_0\): \[f(x_0\pm h) = f(x_0)\pm hf'(x_0)+\frac {h^2}2f''(x_0)\pm \frac {h^3}{6}f'''(\xi _\pm ) .\] Subtracting, the even terms cancel: \(f(x_0+h)-f(x_0-h) = 2hf'(x_0)+\frac {h^3}{6}\bigl (f'''(\xi _+)+f'''(\xi _-)\bigr )\). By the Intermediate Value Theorem the average of \(f'''(\xi _+)\) and \(f'''(\xi _-)\) is \(f'''(\xi _1)\) for some \(\xi _1\), and dividing by \(2h\) gives the formula. The endpoint formula follows in the same way from the expansions of \(f(x_0+h)\) and \(f(x_0+2h)\), choosing the combination that cancels the \(f''\) terms. □
The midpoint formula has half the error constant of the endpoint formula and uses one fewer new value, so it is preferred whenever values on both sides are available; the endpoint formula is for the ends of a table. Using five points gives the five-point midpoint formula \[f'(x_0) = \frac {1}{12h}\bigl [f(x_0-2h)-8f(x_0-h)+8f(x_0+h)-f(x_0+2h)\bigr ]+\frac {h^4}{30}f^{(5)}(\xi ),\] which is \(O(h^4)\).
Example 8.1.2. Values of \(f(x) = xe^x\) are given below. Approximate \(f'(2.0) = 22.167168\) by each formula.
| \(x\) | \(1.8\) | \(1.9\) | \(2.0\) | \(2.1\) | \(2.2\) |
| \(f(x)\) | \(10.889365\) | \(12.703199\) | \(14.778112\) | \(17.148957\) | \(19.855030\) |
Solution.
| formula | approximation | error |
| forward, \(h = 0.1\) | \(\frac {17.148957-14.778112}{0.1} = 23.708450\) | \(1.5\times 10^{0}\) |
| backward, \(h = 0.1\) | \(\frac {14.778112-12.703199}{0.1} = 20.749130\) | \(1.4\times 10^{0}\) |
| three-point endpoint, \(h = 0.1\) | \(22.032310\) | \(1.3\times 10^{-1}\) |
| three-point endpoint, \(h = -0.1\) | \(22.054525\) | \(1.1\times 10^{-1}\) |
| three-point midpoint, \(h = 0.2\) | \(22.414163\) | \(2.5\times 10^{-1}\) |
| three-point midpoint, \(h = 0.1\) | \(22.228790\) | \(6.2\times 10^{-2}\) |
| five-point midpoint, \(h = 0.1\) | \(22.166999\) | \(1.7\times 10^{-4}\) |
Halving \(h\) in the midpoint formula reduced the error by a factor of \(4\), as \(O(h^2)\) predicts, and the five-point formula is several hundred times more accurate than the three-point ones from the same data.
Round-off error
Numerical differentiation is unstable: the formulas divide a small difference of nearly equal function values by the small number \(h\). Suppose each value is computed with an error of at most \(\varepsilon \) and \(|f'''|\leq M\). Then the computed midpoint quotient differs from \(f'(x_0)\) by at most \begin {equation} \label {eq:diffround} e(h) = \frac {\varepsilon }{h}+\frac {h^2}{6}M . \end {equation} The first term grows as \(h\) shrinks; the second shrinks. Setting \(e'(h) = 0\) gives the best step \(h = \left (\frac {3\varepsilon }{M}\right )^{1/3}\): making \(h\) smaller than this makes the result worse.
Example 8.1.3. Approximate \(f'(0.900)\) for \(f(x) = \sin x\) by the three-point midpoint formula, using values of \(\sin x\) rounded to five decimal places.
Solution. The true value is \(\cos 0.9 = 0.62161\).
| \(h\) | \(\sin (0.9+h)\) | \(\sin (0.9-h)\) | approximation | error |
| \(0.1\) | \(0.84147\) | \(0.71736\) | \(0.62055\) | \(-1.06\times 10^{-3}\) |
| \(0.05\) | \(0.81342\) | \(0.75128\) | \(0.62140\) | \(-2.1\times 10^{-4}\) |
| \(0.01\) | \(0.78950\) | \(0.77707\) | \(0.62150\) | \(-1.1\times 10^{-4}\) |
| \(0.005\) | \(0.78643\) | \(0.78021\) | \(0.62200\) | \(3.9\times 10^{-4}\) |
| \(0.001\) | \(0.78395\) | \(0.78270\) | \(0.62500\) | \(3.4\times 10^{-3}\) |
The error is smallest for \(h\) between \(0.01\) and \(0.05\) and then grows again. With \(\varepsilon = 5\times 10^{-6}\) and \(M = \max |{-\cos x}|\leq 1\), the theory gives the optimal step \(h = (3\times 5\times 10^{-6})^{1/3} = 0.025\) (Figure 9, right).
Richardson extrapolation
If an approximation \(N(h)\) of a quantity \(M\) has an error with a known form, approximations for two values of \(h\) can be combined to cancel the leading error term. Suppose \[M = N_1(h)+K_1h^2+K_2h^4+\cdots .\] Replacing \(h\) by \(\frac h2\) gives \(M = N_1(\frac h2)+K_1\frac {h^2}{4}+K_2\frac {h^4}{16}+\cdots \). Multiplying this by \(4\), subtracting the first equation and dividing by \(3\) eliminates the \(h^2\) term: \[M = N_2(h)+O(h^4),\qquad N_2(h) = N_1\!\left (\frac h2\right )+\frac {N_1(\frac h2)-N_1(h)}{3} .\] This is Richardson extrapolation. It can be repeated: if only even powers of \(h\) occur, \(N_{j}(h) = N_{j-1}(\frac h2)+\frac {N_{j-1}(\frac h2)-N_{j-1}(h)}{4^{j-1}-1}\) removes the next term each time.
Example 8.1.4. Combine the two midpoint approximations of Example 8.1.2 by Richardson extrapolation.
Solution. The midpoint formula has an error expansion in even powers of \(h\) (by the symmetry of the Taylor series), and \(N_1(0.2) = 22.414163\), \(N_1(0.1) = 22.228790\). So \[N_2(0.2) = 22.228790+\frac {22.228790-22.414163}{3} = 22.166999,\] with error \(1.7\times 10^{-4}\): the accuracy of the five-point formula, obtained by combining two three-point results. (In fact this combination is the five-point formula.)
Questions on this section
Stuck on something here? Ask below and it stays attached to this topic.