Skip to main content

Numerical Analysis A look at rootfinding, quadrature and approximation using analytic and computing tools

Peter Staab

Section 5.1 Numerical Differentiation

In this section, we examine how to approximate the derivative of a function at some point. There are numerous approximations and we will also estimate the error.
🔗
Recall that the definition of the derivative of a function at a point \(x\) is
\begin{equation*} f'(x) = \lim_{h \rightarrow 0}\frac{f(x+h)-f(x)}{h} \end{equation*}
So for some \(h \ge 0\text{,}\) we can write down an estimate of the derivative with
\begin{align} f'(x_{0})\amp \approx \frac{f(x_{0}+h)-f(x_{0})}{h} \notag\\ \amp = \frac{f(x_{1})-f(x_{0})}{h} \tag{5.1.1} \end{align}
where \(x_{1} =x_{0}+h\text{.}\) This formula is called the forward difference formula. Note also that this formula is
\begin{equation*} f[x_{0},x_{1}] = \frac{f(x_{1})-f(x_{0})}{x_{1}-x_{0}} \end{equation*}
the Newton Divided Difference formula as in (4.3.5). Thus, as we mentioned in Section 4.3 divided difference formulas are like derivatives.
🔗

Subsection 5.1.1 Two-Point Difference Formulas

In this section, we examine the derivative formula as well as errors associated with the derivative. First, we will consider derivative formulas that are based on two points.
🔗
Let
\begin{equation} f(x) = P_{1}(x) + \frac{f''(\xi(x))}{2!}(x-x_{0})(x-x_{1})\tag{5.1.2} \end{equation}
where \(P_{1}(x)\) is the interpolating polynomial that passes through \((x_{0},f(x_{0}))\) and \((x_{1},f(x_{1}))\) and the second term is the error (as seen in Theorem 4.2.6).
🔗
Since we can write \(P_{1}(x)\) as the Lagrange form of the linear interpolating polynomial from (4.2.2) and (4.2.3) as
\begin{equation*} \begin{aligned}P_{1}(x)\amp = f(x_{0}) \frac{x-x_{1}}{x_{0}-x_{1}}+ f(x_{1}) \frac{x-x_{0}}{x_{1}-x_{0}}\\\amp = \frac{f(x_{0})}{-h}(x-x_{1}) + \frac{f(x_{1})}{h}(x-x_{0})\end{aligned} \end{equation*}
and differentiating, the results are
\begin{equation*} \begin{aligned}P_{1}'(x)\amp = \frac{f(x_{1})-f(x_{0})}{h}.\end{aligned} \end{equation*}
🔗
Differentiating (5.1.2), one gets
\begin{align} f'(x)\amp = P_{1}'(x) + D_{x} \biggl(\frac{f''(\xi(x))}{2!}(x-x_{0})(x-x_{1}) \biggr) \notag\\ \amp \qquad \text{Differentiating the second term using the product rule,} \notag\\ \amp = \frac{f(x_{1})-f(x_{0})}{h}+ D_{x} \biggl(\frac{f''(\xi(x))}{2}\biggr) (x-x_{0})(x-x_{1}) \notag\\ \amp \qquad + \frac{f''(\xi(x))}{2}( 2x -x_{0} -x_{1} ) \tag{5.1.3} \end{align}
🔗
Evaluating (5.1.3) at both \(x=x_{0}\) and \(x=x_{1}\text{,}\) one gets:
\begin{equation*} \begin{aligned}f'(x_{0})\amp = P_{1}'(x_{0}) + \frac{h}{2}f''(\xi) = \frac{f(x_{1})-f(x_{0})}{h}+ \frac{h}{2}f''(\xi) \\ f'(x_{1})\amp =P_{1}'(x_{1}) + \frac{h}{2}f''(\xi) = \frac{f(x_{1})-f(x_{0})}{h}+ \frac{h}{2}f''(\xi).\end{aligned} \end{equation*}
The forward difference formula and backwards difference formula for \(f'(x)\) are respectively:
\begin{equation*} \begin{aligned}f'(x_{0})\amp \approx \frac{f(x_{1})-f(x_{0})}{h}\\ f'(x_{1})\amp \approx \frac{f(x_{1})-f(x_{0})}{h}\end{aligned} \end{equation*}
and the error associated with each is \(O(h)\text{.}\) In particular, the error is \((h/2)f''(\xi)\) for some \(\xi \in [x_{0},x_{1}]\text{.}\)
🔗
In the case of the 2-point first derivatives, the forward and backward difference formulas are the same. We will see for 3-point (and higher) formulas that they won’t be.
🔗
🔗

Subsection 5.1.2 \(n\)-point formula for \(f'(x)\)

Now that we have seen the 2-point formulas, we can generalize to \(n\)-point formulas using the following steps. If the function is known at the \(n\) points: \(x_{0}\text{,}\) \(x_{1}\text{,}\) \(\ldots\text{,}\) \(x_{n-1}\text{,}\) then to find \(f'(x)\) at these points, then
🔗
  1. Write the general interpolating polynomial, \(P_{n-1}(x)\text{.}\) through the \(n\)-points \(x_{0},x_{1}, \ldots, x_{n-1}\text{.}\) (Generally the Lagrange form of the polynomial is used.)
    🔗
    🔗
  2. Find \(P'_{n-1}(x)\text{.}\)
    🔗
    🔗
  3. Evaluate \(P'_{n-1}(x_{i})\) for \(i=0,1,\ldots, n-1\text{.}\)
    🔗
    🔗
🔗
Also, generally, the formulas are found such that \(x_{i+1}-x_{i}=h\text{,}\) that is the points are equally-spaced with spacing \(h\text{.}\) If, in addition, the error formula is sought, recall that the polynomial interpolation error \(E(x)\) is given as
\begin{equation*} \begin{aligned}f(x)\amp = P_{n-1}(x) + E_{n-1}(x)\end{aligned} \end{equation*}
where
\begin{equation*} \begin{aligned}E_{n-1}(x)\amp = \frac{f^{(n)}(\xi(x))}{n!}(x-x_{0})(x-x_{1}) \cdots(x-x_{n-1})\end{aligned} \end{equation*}
🔗
Then find \(E'_{n-1}(x_{i})\) for the desired point \(x_{i}\text{.}\)
🔗
🔗

Subsection 5.1.3 3- and 5-point formulas for \(f'(x)\)

We now apply the general results to find 3-point formulas for \(f'(x)\text{.}\) First, to find the 3-point formula we find the interpolated polynomial \(P_{2}(x)\text{:}\)
\begin{equation*} \begin{aligned}P_{2}(x)\amp = f(x_{0}) \frac{(x-x_{1})(x-x_{2})}{(x_{0} - x_{1})(x_{0} -x_{2})}+f(x_{1}) \frac{(x-x_{0})(x-x_{2})}{(x_{1} - x_{0})(x_{1} -x_{2})}, \\\amp \qquad +f(x_{2}) \frac{(x-x_{0})(x-x_{1})}{(x_{2} - x_{0})(x_{2} -x_{1})},\end{aligned} \end{equation*}
and let \(h=x_{1}-x_{0}\) and \(h=x_{2}-x_{1}\text{,}\)
\begin{align} P_{2} (x)\amp = \frac{f(x_{0})}{(-h)(-2h)}(x-x_{1})(x-x_{2}) + \frac{f(x_{1})}{h(-h)}(x-x_{0})(x-x_{2}) + \frac{f(x_{2})}{(2h)(h)}, \notag\\ \amp = \frac{1}{2h^{2}}\bigl(f(x_{0}) (x-x_{1})(x-x_{2}) -2f(x_{1}) (x-x_{0})(x-x_{2}), \notag\\ \amp \qquad + f(x_{2}) (x-x_{0})(x-x_{1}) \bigr), \notag\\ \amp \qquad \quad \text{and differentiating, one gets } \notag\\ P'_{2}(x)\amp = \frac{1}{2h^{2}}\bigl( f(x_{0}) [(x-x_{1}) + (x-x_{2})] - 2 f(x_{1}) [ (x-x_{0}) + (x-x_{2}) ], \notag\\ \amp \qquad \qquad + f(x_{2}) [ (x-x_{0}) + (x-x_{1}) ] \bigr). \tag{5.1.4} \end{align}
🔗
There are three formulas for 3-pt difference formulas when (5.1.4) is evaluated \(x_{0}\text{,}\) \(x_{1}\) and \(x_{2}\text{.}\)
\begin{align} P'_{2}(x_{0})\amp = \frac{1}{2h^{2}}\bigl( -3h f(x_{0}) + 4h f(x_{1}) - h f(x_{2}) \bigr),\notag\\ \amp = \frac{1}{2h}\bigl( -3 f(x_{0}) +4 f(x_{1}) - f(x_{2}) \bigr). \tag{5.1.5}\\ P'_{2}(x_{1})\amp = \frac{1}{2h^{2}}\bigl(-h f(x_{0}) + h f(x_{2}) \bigr), \notag\\ \amp = \frac{1}{2h}\bigl( f(x_{2}) -f(x_{0}) \bigr). \tag{5.1.6}\\ P'_{2}(x_{2})\amp = \frac{1}{2h^{2}}\bigl( h f(x_{0}) - 2 (2h) f(x_{1}) + (3h) f(x_{2}) \bigr), \notag\\ \amp = \frac{1}{2h}\bigl( f(x_{0}) - 4 f(x_{1}) + 3 f(x_{2}) \bigr), \tag{5.1.7} \end{align}
where (5.1.5) is called the 3-point forward difference formula for \(f'(x)\), (5.1.6) is called the 3-point center difference formula for \(f'(x)\) and (5.1.7) is called the 3-point backward difference formula for \(f'(x)\).
🔗
Next, we examine the error formulas for these numerical derivatives. The error formula for the order-2 interpolation is:
\begin{equation*} E_{2}(x) = \frac{f'''(\xi(x))}{3!}(x-x_{0})(x-x_{1})(x-x_{2}) \end{equation*}
and differentiating leads to
\begin{align} E_{2}'(x)\amp = D_{x} \biggl( \frac{f'''(\xi(x))}{6}(x-x_{0})(x-x_{1})(x-x_{2}), \biggr) \notag\\ \amp \qquad \text{And using the product rule results in} \notag\\ \amp = D_{x} \biggl(\frac{f'''(\xi(x))}{6}\biggr) (x-x_{0})(x-x_{1})(x-x_{2}), \notag\\ \amp \qquad + \frac{f'''(\xi(x))}{6}\bigl( (x-x_{1})(x-x_{2}) + (x-x_{0})(x-x_{2}), \notag\\ \amp \qquad + (x-x_{0})(x-x_{1}) \bigr). \tag{5.1.8} \end{align}
🔗
To find the error associated with the forward, center and backward difference formulas, evaluate (5.1.8) at \(x_{0}, x_{1}\) and \(x_{2}\) to get
\begin{equation*} \begin{aligned}E_{2}'(x_{0})\amp = 0 + \frac{f'''(\xi)}{6}(2h^{2}) = \frac{h^{2}}{3}f'''(\xi) \\ E_{2}'(x_{1})\amp = 0 + \frac{f'''(\xi)}{6}(-h^{2}) = - \frac{h^{2}}{6}f'''(\xi) \\ E_{2}'(x_{2})\amp = 0 + \frac{f'''(\xi)}{6}(2h^{2}) = \frac{h^{2}}{3}f'''(\xi)\end{aligned} \end{equation*}
🔗
The main importance of these formulas is that they are all \(O(h^{2})\) methods.
🔗

Example 5.1.1.

Use the 3-point forward difference formula in (5.1.5) above to estimate \(f'(0)\) if \(f(x) = \sin x\) and \(x_{0}=0, x_{1}=\pi/6\) and \(x_{2} =\pi/3\text{.}\) Provide a theoretical error bound and compare it to the actual derivative.
🔗
Solution.
First, the forward difference formula applied to this problem is
\begin{equation*} \begin{aligned}f'(x_{0})\amp \approx \frac{1}{2h}\bigl(- 3 f(0) + 4 f(\pi/6) - f(\pi/3) \bigr) \\\amp = \frac{1}{2\pi/6}\bigl( 0 + 2 - \frac{\sqrt{3}}{2}\bigr) \\\amp \approx 1.083\end{aligned} \end{equation*}
and recall that \(f'(0)=\cos 0 =1\text{.}\)
🔗
The theoretical error bound is
\begin{equation*} \begin{aligned}|E_{2}'(x_{0})|\amp = \frac{h^{2}}{3}|f'''(\xi)| \\ \amp \qquad \text{and since $f(x)=\sin x$, $|f'''(\xi)| \leq 1$}\\ \amp \leq \biggl(\frac{\pi}{6}\biggr)^{2} \frac{1}{3}\approx 0.0913\end{aligned} \end{equation*}
and since the actual error is 0.083, within the theoretical error bound.
🔗
🔗
🔗
The next example uses a CAS to do the algebra and calculus to find the 5-point formula for \(f'(x)\text{.}\)
🔗

Example 5.1.2.

Find the 5-point center difference formula for \(f'(x)\) using a Computer Algebra System:
🔗
Solution.
We will use the Symbolics package of Julia.
 1 
This was last done in March of 2026 using Julia 1.13 and the Symbolics 7.16.0 package. It appears that Symbolics is changing quickly and may be able to do more with newer packages or the syntax may change.
🔗
using Symbolics
🔗
Then set up the needed variables.
🔗
@variables x[0:4], h, f(..), t
🔗
Then we need to define the order-4 Lagrange polynomial through the points \((x_k, f(x_k))\) for \(k=0,1,\ldots, 4\text{.}\) Unfortunately, it looks like it needs to be written out:
🔗
     P4 = f(x[0]) * (t-x[1])*(t-x[2])*(t-x[3])*(t-x[4])/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3])*(x[0]-x[4])) + 
f(x[1])*(t-x[0])*(t-x[2])*(t-x[3])*(t-x[4])/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3])*(x[1]-x[4])) +
f(x[2])*(t-x[0])*(t-x[1])*(t-x[3])*(t-x[4])/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3])*(x[2]-x[4])) +
f(x[3])*(t-x[0])*(t-x[1])*(t-x[2])*(t-x[4])/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2])*(x[3]-x[4])) +
f(x[4])*(t-x[0])*(t-x[1])*(t-x[2])*(t-x[3])/((x[4]-x[0])*(x[4]-x[1])*(x[4]-x[2])*(x[4]-x[3]))
and then differentiate:
🔗
dP4 = Symbolics.derivative(P4,t)
which results in a large mess which we won’t show here. To find the difference formula, we substitute in \(x_2\) for \(t\) because we are looking for the center difference. This can be down with
🔗
substitute(substitute(substitute(substitute(dP4, Dict(t => x[2], x[4] => x[3] + h)), 
  x[3] => x[2] + h), 
    x[2] => x[1] + h), 
      x[1] => x[0] + h)
This results in the following:
\begin{equation*} \frac{2 ~ f\left( x_{3} \right)}{3 ~ h} + \frac{ - f\left( x_{4} \right)}{12 ~ h} + \frac{ - 2 ~ f\left( x_{1} \right)}{3 ~ h} + \frac{f\left( x_{0} \right)}{12 ~ h} \end{equation*}
🔗
This can be written as
\begin{equation*} f'(x_2) \approx \frac{1}{12h}\left(f(x_{0})-8f(x_{1})+8f(x_{3})-f(x_{4}) \right) \end{equation*}
🔗
This is the 5-point difference formula for \(f'(x)\) at \(x=x_{2}\text{,}\) the center point. Note that \(f(x_{2})\) is missing from the formula. This is similar to the center point missing in the 3-point formula for \(f'(x)\) in (5.1.6).
🔗
🔗
🔗
🔗

Subsection 5.1.4 Three-point numerical formulas for \(f''(x)\)

For the 3-point formula for \(f''(x)\text{,}\) we need to differentiate the order-2 interpolated polynomial, \(P_{2}(x)\text{.}\) Since we have already found \(P_{2}'(x)\text{,}\) we then differentiate (5.1.4),
\begin{align} P''_{2}(x)\amp = \frac{1}{2h^{2}}\bigl( 2f(x_{0}) - 4 f(x_{1}) +2 f(x_{2}) \bigr) \notag\\ \amp = \frac{1}{h^{2}}\bigl( f(x_{0}) - 2f(x_{1}) + f(x_{2}) \bigr) \tag{5.1.9} \end{align}
and evaluating this at \(x_{0},x_{1}\) and \(x_{2}\) results in the same formula for each. Thus the 3-point forward difference forward for \(f''(x)\), the 3-point center difference forward for \(f''(x)\), and the 3-point backward difference forward for \(f''(x)\) are each the same as given in (5.1.9).
🔗
To determine the error in this formula, differentiate (5.1.8).
\begin{equation*} \begin{aligned}E_{2}''(x)\amp = D^{2}_{x} \biggl(\frac{f'''(\xi(x))}{6}\biggr) (x-x_{0})(x-x_{1})(x-x_{2}) \notag \\\amp \qquad + D_{x} \biggl( \frac{f'''(\xi(x))}{6}\biggr) \bigl( (x-x_{1})(x-x_{2}) + (x-x_{0})(x-x_{2}) \\\amp \qquad \qquad + (x-x_{0})(x-x_{1}) \bigr) \\\amp \qquad + \frac{f'''(\xi(x))}{6}\bigl( 2(x-x_{0}) + 2(x-x_{1}) + 2(x-x_{2}) \bigr) \\\end{aligned} \end{equation*}
and evaluating this at \(x_{0}, x_{1}\) and \(x_{2}\) results in
\begin{equation*} \begin{aligned}E_{2}''(x_{0})\amp = D_{x} \biggl( \frac{f'''(\xi(x))}{6}\biggr) (2h^{2}) + 6h \frac{f'''(\xi(x))}{6}\\ E_{2}''(x_{1})\amp = D_{x} \biggl( \frac{f'''(\xi(x))}{6}\biggr) (-h^{2}) \\ E_{2}''(x_{2})\amp = D_{x} \biggl( \frac{f'''(\xi(x))}{6}\biggr) (2h^{2}) + 6h \frac{f'''(\xi(x))}{6}\\\end{aligned} \end{equation*}
🔗
Thus the forward and backward difference formulas are \(O(h)\) methods and the center difference formula is \(O(h^{2})\text{.}\)
🔗
🔗

Subsection 5.1.5 Higher Order Numerical derivatives

The examples in this section show how to derive specific formulas for \(f'\) and \(f''\text{.}\) Following similar techniques, one can derive forward, center and backward difference formulas for any order derivative using the same basic steps:
🔗
  1. Find an interpolating polynomial, \(P_{n}(x)\) that passes through the desired points. (Note: often these are equally spaced in that \(x_{i+1}-x_{i}=h\text{,}\) however this is not necessary.) The order of the polynomial must be at least one less than the number of points in the formula (That is, a 7-point formula must have a order of at least 6).
    🔗
    🔗
  2. Find the desired derivative of \(P_{n}(x)\text{.}\)
    🔗
    🔗
  3. Substitute in the value of \(x\) where one wants the derivative.
    🔗
    🔗
🔗
🔗
🔗
PrevTopNext