from numpy import *
from scipy import linalg
from scipy.linalg import norm
from matplotlib.pyplot import *
from prettytable import PrettyTable
from timeit import default_timer as timer
import sys
sys.path.append('fncbook/')
import fncbook as FNC

# This (optional) block is for improving the display of plots.
# from IPython.display import set_matplotlib_formats
# set_matplotlib_formats("svg","pdf")
# %config InlineBackend.figure_format = 'svg'
rcParams["figure.figsize"] = [7, 4]
rcParams["lines.linewidth"] = 2
rcParams["lines.markersize"] = 4
rcParams['animation.html'] = "jshtml"  # or try "html5"

5.5.Convergence of finite differences¶

All the finite-difference formulas in the previous section based on equally spaced nodes converge as the node spacing $h$ decreases to zero. However, note that to discretize a function over an interval $[a,b]$ , we use $h=(b-a)/n$ , which implies $n=(b-a)/h=O(h^{-1})$ . As $h\to 0$ , the total number of nodes needed grows without bound. So we would like to make $h$ as large as possible while still achieving some acceptable accuracy.

Although we are measuring the truncation error only at $x=0$ , it could be defined for other $x$ as well. The definition adjusts naturally to use $f''(0)$ for difference formulas targeting the second derivative.

All the finite-difference formulas given in Finite differences are convergent.

5.5.1Order of accuracy¶

Of major interest is the rate at which $\tau_f\to 0$ in a convergent formula.

Hence, the forward-difference formula in Example 5.5.1 has order of accuracy equal to 1; i.e., it is first-order accurate. All else being equal, a higher order of accuracy is preferred, since $O(h^m)$ vanishes more quickly for larger values of $m$ . As a rule, including more function values in a finite-difference formula (i.e., increasing the number of weights in (5.4.1)) increases the order of accuracy, as can be seen in Table 5.4.1 and Table 5.4.2.

Order of accuracy is calculated by expanding $\tau_f$ in a Taylor series about $h=0$ and ignoring all but the leading term.^[1]

Example 5.5.2

We compute the truncation error of the centered difference formula (5.4.5):

\begin{split} \tau_f(h) &= f'(0) - \frac{ f(h)-f(-h)}{2h}\\ &= f'(0) - (2h)^{-1} \left[ \bigl( f(0) + h f'(0) + \tfrac{1}{2}h^2f''(0)+ \tfrac{1}{6}h^3f'''(0)+ O(h^4) \bigr) \right.\\ &\qquad - \left. \bigl( f(0) - h f'(0) + \tfrac{1}{2}h^2f''(0) - \tfrac{1}{6}h^3f'''(0)+O(h^4) \bigr) \right] \\ &= -(2h)^{-1} \left[ \tfrac{1}{3}h^3f'''(0) + O(h^4) \right] = O(h^2). \end{split}

(5.5.3)

Thus, this method has order of accuracy equal to 2.

Example 5.5.3 (Convergence of finite differences)

Let’s observe the convergence of the formulas in Example 5.5.1 and Example 5.5.2, applied to the function $\sin(e^{x+1})$ at $x=0$ .

f = lambda x: sin(exp(x + 1))
exact_value = exp(1) * cos(exp(1))

We’ll compute the formulas in parallel for a sequence of $h$ values.

h_ = array([5 / 10**(n+1) for n in range(6)])
FD = zeros((len(h_), 2))
for (i, h) in enumerate(h_):
    FD[i, 0] = (f(h) - f(0)) / h 
    FD[i, 1] = (f(h) - f(-h)) / (2*h)
results = PrettyTable()
results.add_column("h", h_)
results.add_column("FD1", FD[:, 0])
results.add_column("FD2", FD[:, 1])
print(results)

+--------+---------------------+---------------------+
|   h    |         FD1         |         FD2         |
+--------+---------------------+---------------------+
|  0.5   |  -2.768575766550465 | -1.9704719803862911 |
|  0.05  | -2.6127952856136947 |  -2.475520256824717 |
| 0.005  | -2.4921052189334048 | -2.4783216951179465 |
| 0.0005 | -2.4797278609705042 | -2.4783494526022243 |
| 5e-05  | -2.4784875710526233 |  -2.478349730152263 |
| 5e-06  |  -2.478363517077753 |  -2.478349732970564 |
+--------+---------------------+---------------------+

All that’s easy to see from this table is that FD2 appears to converge to the same result as FD1, but more rapidly. A table of errors is more informative.

errors = FD - exact_value
results = PrettyTable()
results.add_column("h", h_)
results.add_column("error in FD1", errors[:, 0])
results.add_column("error in FD2", errors[:, 1])
print(results)

+--------+-------------------------+------------------------+
|   h    |       error in FD1      |      error in FD2      |
+--------+-------------------------+------------------------+
|  0.5   |   -0.29022603359523025  |   0.5078777525689437   |
|  0.05  |   -0.1344455526584598   | 0.0028294761305178717  |
| 0.005  |   -0.01375548597816989  | 2.8037837288330536e-05 |
| 0.0005 |  -0.0013781280152693753 | 2.8035301058437767e-07 |
| 5e-05  |  -0.0001378380973884319 | 2.8029720766653554e-09 |
| 5e-06  | -1.3784122518067932e-05 | -1.532907134560446e-11 |
+--------+-------------------------+------------------------+

In each row, $h$ is decreased by a factor of 10, so that the error is reduced by a factor of 10 in the first-order method and 100 in the second-order method.

A graphical comparison can be useful as well. On a log-log scale, the error should (as $h\to 0$ ) be a straight line whose slope is the order of accuracy. However, it’s conventional in convergence plots to show $h$ decreasing from left to right, which negates the slopes.

plot(h_, abs(errors), "o-", label=["FD1", "FD2"])
gca().invert_xaxis()
# Add lines for perfect 1st and 2nd order.
loglog(h_, h_, "--", label="$O(h)$")
loglog(h_, h_**2, "--", label="$O(h^2)$")
xlabel("$h$")
ylabel("error")
legend();

5.5.2Stability¶

The truncation error $\tau_f(h)$ of a finite-difference formula is dominated by a leading term $O(h^m)$ for an integer $m$ . This error decreases as $h\to 0$ . However, we have not yet accounted for the effects of roundoff error. To keep matters as simple as possible, let’s consider the forward difference

\delta(h) = \frac{f(x+h)-f(x)}{h}.

(5.5.4)

As $h\to 0$ , the numerator approaches zero even though the values $f(x+h)$ and $f(x)$ are not necessarily near zero. This is the recipe for subtractive cancellation error! In fact, finite-difference formulas are inherently ill-conditioned as $h\to 0$ . To be precise, recall that the condition number for the problem of computing $f(x+h)-f(x)$ is

\kappa(h) = \frac{ \max\{\,|f(x+h)|,|f(x)|\,\} }{ |f(x+h)-f(x) | },

(5.5.5)

implying a relative error of size $\kappa(h) \macheps$ in its computation. Hence, the numerical value we actually compute for $\delta$ is

\begin{split} \tilde{\delta}(h) &= \frac{f(x+h)-f(x)}{h}\, (1+\kappa(h)\macheps) \\ &= \delta(h) + \frac{ \max\{\,|f(x+h)|,|f(x)|\,\} }{ |f(x+h)-f(x) | }\cdot \frac{f(x+h)-f(x)}{h} \cdot \macheps. \end{split}

(5.5.6)

Hence, as $h\to 0$ ,

\bigl| \tilde{\delta}(h) - \delta(h) \bigr| = \frac{ \max\{\,|f(x+h)|,|f(x)|\,\} }{ h}\,\macheps \sim |f(x)|\, \macheps\cdot h^{-1}.

(5.5.7)

Combining the truncation error and the roundoff error leads to

\bigl| f'(x) - \tilde{\delta}(h) \bigr| \le \bigl| \tau_f(h) \bigr| + \bigl|f(x) \bigr|\, \macheps \, h^{-1}.

(5.5.8)

Equation (5.5.8) indicates that while the truncation error $\tau$ vanishes as $h$ decreases, the roundoff error actually increases thanks to the subtractive cancellation. At some value of $h$ the two error contributions will be of roughly equal size. This occurs when

\bigl|f(x)\bigr|\, \macheps\, h^{-1} \approx C h, \quad \text{or} \quad h \approx K \sqrt{\rule[0.05em]{0mm}{0.4em}\macheps},

(5.5.9)

for a constant $K$ that depends on $x$ and $f$ , but not $h$ . In summary, for a first-order finite-difference method, the optimum spacing between nodes is proportional to $\macheps^{\,\,1/2}$ . (This observation explains the choice of δ in Function 4.6.1.)

For a method of truncation order $m$ , the details of the subtractive cancellation are a bit different, but the conclusion generalizes.

A different statement of the conclusion is that for a first-order formula, at most we can expect accuracy in only about half of the available machine digits. As $m$ increases, we get ever closer to using the full accuracy available. Higher-order finite-difference methods are both more efficient and less vulnerable to roundoff than low-order methods.

Example 5.5.4 (Roundoff error in finite differences)

Let $f(x)=e^{-1.3x}$ . We apply finite-difference formulas of first, second, and fourth order to estimate $f'(0)=-1.3$ .

f = lambda x: exp(-1.3 * x)
exact = -1.3

h_ = array([1 / 10**(n+1) for n in range(12)])
FD = zeros((len(h_), 3))
for (i, h) in enumerate(h_):
    nodes = h * linspace(-2, 2, 5)
    vals = f(nodes)
    FD[i, 0] = dot(array([0, 0, -1, 1, 0]) / h, vals)
    FD[i, 1] = dot(array([0, -1/2, 0, 1/2, 0]) / h, vals)
    FD[i, 2] = dot(array([1/12, -2/3, 0, 2/3, -1/12]) / h, vals)

results = PrettyTable()
results.add_column("h", h_)
results.add_column("FD1", FD[:, 0])
results.add_column("FD2", FD[:, 1])
results.add_column("FD4", FD[:, 2])
print(results)

+--------+---------------------+---------------------+---------------------+
|   h    |         FD1         |         FD2         |         FD4         |
+--------+---------------------+---------------------+---------------------+
|  0.1   |  -1.219045690794387 | -1.3036647620203026 | -1.2999875986418985 |
|  0.01  | -1.2915864979712421 | -1.3000366169760815 | -1.2999999987623323 |
| 0.001  |  -1.299155366047744 | -1.3000003661667074 | -1.2999999999998653 |
| 0.0001 | -1.2999155036619303 |  -1.300000003662075 | -1.2999999999996963 |
| 1e-05  | -1.2999915500411239 | -1.3000000000396366 | -1.3000000000099614 |
| 1e-06  | -1.2999991549911272 | -1.2999999999900496 |  -1.300000000024955 |
| 1e-07  |  -1.299999915493899 | -1.3000000000289944 | -1.3000000004946557 |
| 1e-08  | -1.2999999965401798 | -1.3000000042305544 | -1.3000000107468235 |
| 1e-09  | -1.2999999965401796 | -1.3000000340328768 | -1.3000000442867095 |
| 1e-10  | -1.3000001075624823 | -1.2999996723115146 | -1.2999999661783477 |
| 1e-11  | -1.3000045484545808 | -1.3000038001061966 | -1.3000094495401033 |
| 1e-12  | -1.2999601395335958 |  -1.300004483829298 | -1.3000696164874246 |
+--------+---------------------+---------------------+---------------------+

They all seem to be converging to -1.3. The convergence plot reveals some interesting structure to the errors, though.

loglog(h_, abs(FD[:, 0] + 1.3), "-o", label="FD1")
loglog(h_, abs(FD[:, 1] + 1.3), "-o", label="FD2")
loglog(h_, abs(FD[:, 2] + 1.3), "-o", label="FD4")
gca().invert_xaxis()
plot(h_, 0.1 * 2 ** (-52) / h_, "--", color="k", label="$O(h^{-1})$")
xlabel("$h$")
ylabel("total error")
title("FD error with roundoff")
legend();

Again the graph is made so that $h$ decreases from left to right. The errors are dominated at first by truncation error, which decreases most rapidly for the fourth-order formula. However, increasing roundoff error eventually equals and then dominates the truncation error as $h$ continues to decrease. As the order of accuracy increases, the crossover point moves to the left (greater efficiency) and down (greater accuracy).

5.5.3Exercises¶

Exercise 5.5.6

✍ A different way to derive finite-difference formulas is the method of undetermined coefficients. Starting from (5.4.1),

f'(x) \approx \frac{1}{h}\sum_{k=-p}^q a_k f(x+kh),

let each $f(x+k h)$ be expanded in a series around $h=0$ . When the coefficients of powers of $h$ are collected, one obtains

\frac{1}{h} \sum_{k=-p}^q a_k f(x+kh) = \frac{b_0}{h} + b_1 f'(x) + b_2 f''(x)h + \cdots,

where

b_i = \sum_{k=-p}^q k^i a_k.

In order to make the result as close as possible to $f'(x)$ , we impose the conditions

b_0 = 0,\, b_1=1,\, b_2=0,\, b_3=0,\,\ldots,\,b_{p+q}=0.

This provides a system of linear equations for the weights.

(a) For $p=q=2$ , write out the system of equations for $a_{-2}$ , $a_{-1}$ , $a_0$ , $a_1$ , $a_2$ .

(b) Verify that the coefficients from the appropriate row of Table 5.4.1 satisfy the equations you wrote down in part (a).

(c) Derive the finite-difference formula for $p=1$ , $q=2$ using the method of undetermined coefficients.

Footnotes¶

The term truncation error is derived from the idea that the finite-difference formula, being finite, has to truncate the series representation and thus cannot be exactly correct for all functions.
↩