Numerical differentiation: Difference formulas for the first derivative
General difference formulas
The difference formulas presented earlier have been introduced via Taylor's theorem. In some applications it is necessary to use difference formulas with a higher accuracy (think of rounding errors) or with grid points that are not equidistantly distributed. Here we will present a method to derive new difference formulas. The 3-point and 5-point central difference formulas are special cases of derivatives computed through a Savitsky-Golay filter.
The general approach for \(n+1\) different grid points \(x_0, x_1, \ldots x_n\), which differ from each other by a multiple of the step size \(h\), is to find coefficients \(c_0, c_1, \ldots, c_n\) such that the expression \[Q(h)\stackrel{\mathrm{def}}{=}\sum_{i=0}^n c_if(x_i)\] approximates the derivative \(f'(x_0)\) with the largest possible order of the truncation error. We find the coefficients \(c_0, c_1, \ldots, c_n\) by developing each \(f(x_i)\) into a Taylor polynomial about \(x_0\) of sufficiently high degree and then by constructing and solving equations for the unknown \(c_i\) 's.
The description above is rather cryptic, but examples make this systematic approach clear.
The 3-point central difference formula We consider 3 grid points \(x_0-h\), \(x_0\) and \(x_0+h\). We define the formula \[Q(h)=A\cdot f(x_0-h)+B\cdot f(x_0)+C\cdot f(x_0+h)\] and determine the Taylor approximation of \(Q(h)\). We know that \[\begin{aligned} f(x_0-h) &= f(x_0)-f'(x_0)h+\tfrac{1}{2}f''(x_0)h^2+O(h^3)\\ f(x_0+h) &= f(x_0)+f'(x_0)h+\tfrac{1}{2}f''(x_0)h^2+O(h^3)\end{aligned}\] and therefore \[\begin{aligned}Q(h) &= A\cdot ( f(x_0)-f'(x_0)h+\tfrac{1}{2}f''(x_0)h^2)+B\cdot f(x_0)\\ &\phantom{=}{}+C\cdot (f(x_0)+f'(x_0)h+\tfrac{1}{2}f''(x_0)h^2) + O(h^3)\\ \\ &= (A+B+C)f(x_0)+(C-A)f'(x_0)h+\tfrac{1}{2}(A+C)f''(x_0)h^2+O(h^3)\end{aligned}\] The equations that \(A, B, C\) must satisfy in order for \(Q(h)\) to maximise the order of the error are \[A+B+C = 0,\quad (C-A)h = 1,\quad A+C = 0\] These equations are easy to solve: \[A=-\frac{1}{2h}, B=0, C=\frac{1}{2h}\] So we have found \[Q(h)=\frac{f(x_0+h)-f(x_0-h)}{2h}\] and \[f'(x_0)=\frac{f(x_0+h)-f(x_0-h)}{2h}+O(h^2)\] This is equal to the central difference formula.
One-sided 3-point difference formulas The left-sided 3-point difference formula: \[f'(x_0)\approx\frac{-3f(x_0)+4f(x_0+h)-f(x_0+2h)}{2h}\] The right-sided 3-point difference formula: \[f'(x_0)\approx\frac{3f(x_0)-4f(x_0-h)+f(x_0-2h)}{2h}\] One-sided 3-point difference formulas are useful for approximating derivatives at the edges of a finite discrete signal.
We take 3 grid points \(x_0\), \(x_0+h\) and \(x_0+2h\). We define the formula \[Q(h)=A\cdot f(x_0)+B\cdot f(x_0+h)+C\cdot f(x_0+2h)\] and determine the Taylor approximation of \(Q(h)\). We know that \[\begin{aligned} f(x_0+h) &= f(x_0)+f'(x_0)h+\tfrac{1}{2}f''(x_0)h^2+O(h^3)\\ f(x_0+2h) &= f(x_0)+2f'(x_0)h+\tfrac{1}{2}f''(x_0)(2h)^2+O(h^3)\end{aligned}\] and therefore \[\begin{aligned}Q(h) &=A\cdot f(x_0) + B\cdot ( f(x_0)+f'(x_0)h+\tfrac{1}{2}f''(x_0)h^2)\\ &\phantom{=}{}+C\cdot (f(x_0)+2f'(x_0)h+\tfrac{4}{2}f''(x_0)h^2) + O(h^3)\\ \\ &= (A+B+C)f(x_0)+(B+2C)f'(x_0)h+\tfrac{1}{2}(B+4C)f''(x_0)h^2+O(h^3)\end{aligned}\] The equations that \(A, B, C\) must satisfy in order for \(Q(h)\) to maximise the order of the error are \[A+B+C = 0,\quad (B+2C)h = 1,\quad B+4C = 0\] These equations are easy to solve: \[A=-\frac{3}{2h}, B=\frac{2}{h}, C=-\frac{1}{2h}\] So we have found \[Q(h)=\frac{-3f(x_0)+4f(x_0+h)-f(x_0+2h)}{2h}\] and \[f'(x_0)=\frac{-3f(x_0)+4f(x_0+h)-f(x_0+2h)}{2h}+O(h^2)\] This is a left-sided 3-point difference formula.
Similarly, you can also find the right-sided 3-point difference formula: \[f'(x_0)=\frac{3f(x_0)-4f(x_0-h)+f(x_0-2h)}{2h}+O(h^2)\]
The 5-point central difference formula \[f'(x_0)\approx\frac{f(x_0-2h)-8f(x_0-h)+8f(x_0+h)-f(x_0+2h)}{12h}\] The truncation error is generally smaller than with the 3-point central difference formula.
We consider 5 grid points \(x_0-2h\), \(x_0-h\), \(x_0\), \(x_0+h\) and \(x_0+2h\). We define the formula \[\begin{aligned}Q(h)&=A\cdot f(x_0-2h)+B\cdot f(x_0-h)+C\cdot f(x_0)\\ &\phantom{=}{}+ D\cdot f(x_0+h)+E\cdot f(x_0+2h)\end{aligned} \] and determine the Taylor approximation of \(Q(h)\). We know that \[\begin{aligned} f(x_0-2h) &= f(x_0)-2f'(x_0)h+2f''(x_0)h^2-\tfrac{4}{3}f'''(x_0)h^3+\tfrac{2}{3}f''''(x_0)h^4+O(h^5)\\ f(x_0-h) &= f(x_0)-f'(x_0)h+\tfrac{1}{2}f''(x_0)h^2-\tfrac{1}{6}f'''(x_0)h^3+\tfrac{1}{24}f''''(x_0)h^4+O(h^5)\\ f(x_0+h) &= f(x_0)+f'(x_0)h+\tfrac{1}{2}f''(x_0)h^2+\tfrac{1}{6}f'''(x_0)h^3+\tfrac{1}{24}f''''(x_0)h^4+O(h^5)\\ f(x_0+2h) &= f(x_0)+2f'(x_0)h+2f''(x_0)h^2+\tfrac{4}{3}f'''(x_0)h^3+\tfrac{2}{3}f''''(x_0)h^4+O(h^5)\end{aligned}\] and therefore \[\begin{aligned}Q(h) &= (A+B+C+D+E)f(x_0)+(-2A-B+D+2E)f'(x_0)h\\&\phantom{=}{}+\tfrac{1}{2}(4A+B+D+4E)f''(x_0)h^2+\tfrac{1}{6}(-8A-B+D+8E)f'''(x_0)h^3\\&\phantom{=}{}+\tfrac{1}{24}(16A+B+D+16E)f''''(x_0)h^4+O(h^5)\end{aligned}\] The equations that \(A, B, C, D, E\) must satisfy in order for \(Q(h)\) to maximise the order of the error are \[\begin{aligned}A+B+C +D+\;\;E&= 0\\ -2A-B+\phantom{+C}D+\;2E &= \frac{1}{h}\\ 4A+B+\phantom{+C}D+\;4E &= 0\\ -8A-B+\phantom{+C}D+\;8E &= 0\\ 16A+B+\phantom{+C}D+16E &= 0\end{aligned}\] These equations are relatively easy to solve: \[A=\frac{1}{12h}, B=-\frac{2}{3h},C=0, D=\frac{2}{3h}, E=-\frac{1}{12h}\] So we have found \[Q(h)=\frac{f(x_0-2h)-8f(x_0-h)+8f(x_0+h)-f(x_0+2h)}{12h}\] and \[f'(x_0)\approx\frac{f(x_0-2h)-8f(x_0-h)+8f(x_0+h)-f(x_0+2h)}{12h}\] This is the 5-point central difference formula. It can be verified that the truncation error has order \(h^4\).