7. Numerical Differentiation and Integration
7.1. Numerical differentiation
In this section we present two methods to derive numerical approximation formulas for the derivative, and we derive some basic approximation formulas.
The derivative of a function is defined by the limit f′(x0) = lim
h→0
f(x0+h)−f(x0)
h .
Therefore, if |h| is small, then the difference quotient f(x0+h)−f(xh 0) is close to the value of the derivative. But we need more: we need to know the truncation error of the approximation. Next we derive this formula in two different ways, and we will derive the formula of the truncation error too.
Supposef ∈C3[a, b] and x0 ∈(a, b). The idea of the first method is the following: We approximate the function f in a neighbourhood of x0 by a Lagrange polynomial Ln(x).
We use L′n(x0) as an approximation of f′(x0). We will call this method as Lagrange’s method. Consider a simple case: let n = 1, x1 =x0 +h ∈(a, b) (and x0 ̸= x1), consider the first-order Lagrange polynomial interpolation of f corresponding to the mesh points x0 and x1:
f(x) =L1(x) +E1(x)
= f(x0)(x−x0−h)
−h +f(x0+h)(x−x0)
h +f′′(ξ(x))
2 (x−x0)(x−x0−h).
Taking the derivative of both sides we get f′(x) = f(x0+h)−f(x0)
h + f′′(ξ(x))
2 (2(x−x0) +h) + d
dx (︂
f′′(ξ(x)))︂(x−x0)(x−x0−h)
2 . (7.1)
Theorem 6.8 yields that the function f′′(ξ(x)) is differentiable forx̸=x0, x0+h, but the derivative cannot be computed explicitly. On the other hand, taking the limit x→x0 in (7.1) we get
f′(x0) = f(x0+h)−f(x0)
h −h
2f′′(ξ), (7.2)
where ξ ∈ ⟨x0, x0+h⟩. Therefore, if we use the approximation formula f′(x0)≈ f(x0+h)−f(x0)
h , (7.3)
the truncation error of the approximation has the form −h2f′′(ξ). Formula (7.3) is called first-order forward difference formula ifh >0, andfirst-order backward difference formula if h < 0. In these formulas the mesh point x0 +h is located right and left to x0, in the respective cases. Formula (7.2) shows that approximation (7.3) is first-order inh. Formula (7.3) is also called two-point difference formula, since it uses two mesh points.
The same formula can be derived (under weaker conditions) in the following way: Let f ∈C2[a, b], and consider the first-order Taylor expansion of f aroundx0:
f(x) =f(x0) +f′(x0)(x−x0) + f′′(ξ(x))
2 (x−x0)2. Substitution x=x0+h gives
f(x0+h) = f(x0) +f′(x0)h+ f′′(ξ) 2 h2, hence
f′(x0) = f(x0+h)−f(x0)
h −h
2f′′(ξ), where ξ =ξ(x0+h).
Example 7.1. Consider the function f(x) = ex2+x. We have f′(x) = ex2+x(2x+ 1), so f′(0) = 1. We compute an approximate value off′(0) using the first-order forward (h >0) and backward (h <0) difference formula, i.e., formula (7.3). In Table 7.1 we printed the approximate values and their errors for different values of h. The numerical results show that if the step size h decreases by one order of magnitude, then the corresponding error also decreases by one order
of magnitude.
□
Table 7.1: First-order difference formula, f(x) = ex2+x, x0 = 0
|h| forward difference error backward difference error 0.100 1.1627807 1.6278e-01 0.8606881 1.3931e-01 0.010 1.0151177 1.5118e-02 0.9851156 1.4884e-02 0.001 1.0015012 1.5012e-03 0.9985012 1.4988e-03
The previous two methods are appropriate to derive higher order, so more precise for-mulas. Suppose f ∈Cn+1, and consider an approximation off by a Lagrange polynomial of degree n:
f(x) =
n
∑︂
k=0
f(xk)lk(x) + f(n+1)(ξ(x))
(n+ 1)! (x−x0)(x−x1)· · ·(x−xn), (7.4)
7.1. Numerical differentiation 139 where lk(x) are the Lagrange basis polynomials of degree n defined by (6.2). Differenti-ating (7.4) and using substitution x=xi we get
f′(xi) =
which is called n+ 1-point difference formula to approximate f′(xi). We apply relation (7.5) for equidistant mesh points, so we assume xj = x0 +jh, where h > 0. It can be shown that the error term in (7.5) is of nth-order in h, and then the resulting formula will also be called difference formula of order n.
Consider the case when n= 2, i.e., we study three-point formulas. Consider the mesh points x0, x0+h, x0+ 2h. Then Relation (7.9) is called three-point midpoint formula orsecond-order central difference formula. (It is also calledcentered difference.) Formula (7.6) is calledthree-point endpoint
formula. It is also called second-order forward difference formula if h > 0, and second-order backward difference formula if h <0.
Example 7.2. We approximate the derivative of the function f(x) = ex2+x at x = 0 with second-order difference formulas (formulas (7.6) and (7.9)). The results can be seen in Table 7.2 for different values of h. The numerical results demonstrate that the truncation error of the
formulas is second-order in h.
□
Table 7.2: Second-order difference formulas, f(x) = ex2+x, x0 = 0
|h| forward error backward error central error 0.100 0.9693157 3.0684e-02 0.9820952 1.7905e-02 1.0117344 1.1734e-02 0.010 0.9997603 2.3968e-04 0.9997728 2.2718e-04 1.0001167 1.1667e-04 0.001 0.9999977 2.3396e-06 0.9999977 2.3271e-06 1.0000012 1.1667e-06
Without proofs we present 5-point central and one-sided formulas, i.e., fourth-order difference formulas:
f′(x0) = 1 12h
(︃
−25f(x0) + 48f(x0+h)−36f(x0+ 2h) + 16f(x0+ 3h)
− 3f(x0+ 4h) )︃
+ h4
5 f(5)(ξ0), (7.10)
f′(x0) = 1 12h
(︃
f(x0−2h)−8f(x0−h) + 8f(x0+h)−f(x0+ 2h) )︃
+ h4
30f(5)(ξ1). (7.11)
Formula (7.10) is one-sided, and (7.11) is central difference.
Example 7.3. We apply formulas (7.10) and (7.11) to approximate the first derivative of f(x) =ex2+x atx= 0. Table 7.3 shows the numerical results.
□
Table 7.3: Fourth-order difference formulas, f(x) = ex2+x, x0 = 0
|h| forward error backward error central error 0.100 0.9967110 3.2890e-03 0.9991793 8.2070e-04 0.9997248 2.7523e-04 0.010 0.9999998 1.7345e-07 0.9999998 1.5136e-07 1.0000000 2.7005e-08 0.001 1.0000000 1.6311e-11 1.0000000 1.6090e-11 1.0000000 2.7000e-12
Next we use the Taylor’s method to derive approximation formulas for higher order derivatives. Let f ∈ C4, and consider the third-order Taylor polynomial expansion of f at x0 with the fourth-order error term:
f(x) = f(x0) +f′(x0)(x−x0) + f′′(x0)
2 (x−x0)2 +f′′′(x0)
6 (x−x0)3+f(4)(ξ)
24 (x−x0)4.
7.1. Numerical differentiation 141 If we substitute x=x0−h and x=x0+h into this relation, we get
f(x0−h) = f(x0)−f′(x0)h+f′′(x0)
2 h2− f′′′(x0)
6 h3+ f(4)(ξ1) 24 h4 and
f(x0+h) = f(x0) +f′(x0)h+f′′(x0)
2 h2+ f′′′(x0)
6 h3+ f(4)(ξ2) 24 h4. Adding the two equations we get
f(x0−h) +f(x0+h) = 2f(x0) +f′′(x0)h2+ f(4)(ξ1) +f(4)(ξ2) 24 h4, which yields
f′′(x0) = f(x0−h)−2f(x0) +f(x0+h)
h2 +f(4)(ξ1) +f(4)(ξ2) 24 h4. Therefore, the approximation formula
f′′(x0)≈ f(x0−h)−2f(x0) +f(x0+h)
h2 (7.12)
has an error of orderh2. We can rewrite the error term f(4)(ξ1)+f24 (4)(ξ2)h4 in a simpler form.
We have by our assumptions that f(4) is continuous, therefore, Theorem 2.2 yields that there exists a point ξ in between ξ1 and ξ2 such that
f(4)(ξ) = f(4)(ξ1) +f(4)(ξ2)
2 .
Hence
f′′(x0) = f(x0−h)−2f(x0) +f(x0+h)
h2 +f(4)(ξ)
12 h2. (7.13)
Example 7.4. We computed the approximation of the second-order derivative off(x) =ex2+x at x = 0 using formula (7.12) and different step sizes. The numerical results can be seen in Table 7.4. Note that the exact derivative value is f′′(0) = 3.
□
Table 7.4: Approximation of the second-order derivative, f(x) = ex2+x, x0 = 0 h approximation error
0.100 3.0209256 2.0926e-02 0.010 3.0002083 2.0834e-04 0.001 3.0000021 2.0833e-06
The numerical differentiation is an unstable problem. To illustrate it we consider a function f(x) and its perturbation of the form
g(x) =f(x) + 1
nsin(n2x).
If we compute an approximation of g′ instead of f′ using any difference formula obtained above, then there is a small change in the function values used in the difference formula if n is large. But the difference between the exact value of the derivatives is large, since g′(x) =f′(x) +ncos(n2x).
Next we investigate the effect of the rounding in numerical differentiation. Consider the simplest difference formula, the first-order difference (7.2). Suppose that here, instead of the exact function values f(x0) and f(x0+h), we use their approximate values f0 and f1, where
f(x0) =f0+e0 and f(x0+h) =f1+e1. Then
f′(x0)≈ f1−f0
h ,
and the resulting error is f′(x0)− f1−f0
h =f′(x0)− f(x0 +h)−f(x0)
h + f(x0+h)−f(x0)
h −f1−f0 h
=−h
2f′′(ξ) + e1 −e0
h . (7.14)
Relation (7.14) shows that the error consists of two parts: the truncation error and the rounding error. If the step size h is small, then the rounding error will be small, but the rounding error can go to ∞as h→0.
Example 7.5. Consider the function f(x) =ex. We compute the approximation off′(1) =e using first-order forward difference formula. In order to enlarge the effect of the rounding, we used 6- and 4-digit arithmetic in the computation. We can see in Table 7.5 that in case of the 4-digit arithmetic, when we decreased the step size to 0.001 from 0.01, the error of the approximation increased. The reason is, clearly, the increase of the rounding error, since here we subtracted two numbers which are close to each other, and also divided by a small number.
□
Table 7.5: Effect of rounding in first-order forward difference, f(x) = ex,x0 = 1 6-digit arithmetic 4-digit arithmetic
h approximation error approximation error 0.100 2.8589000 1.4062e-01 2.8600000 1.4172e-01 0.010 2.7320000 1.3718e-02 2.8000000 8.1718e-02 0.001 2.7200000 1.7182e-03 3.0000000 2.8172e-01
The formulas derived in this section can be applied to approximate partial derivatives.
We list some formulas next.