Chapter 5 - Fundamentals of Numerical Computation

Python implementations

Functions¶

Hat function

hatfun.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
def hatfun(x, t, k):
    """
    hatfun(x,t,k)

    Evaluate a piecewise linear "hat" function at `x`, where `t` is a vector of
    n+1 interpolation nodes and `k` is an integer in 0:n giving the index of the node
    where the hat function equals one.
    """
    n = len(t) - 1

    # Return correct node given mathematical index k, including fictitious choices.
    def node(k):
        if k < 0:
            return 2 * t[0] - t[1]
        elif k > n:
            return 2 * t[n] - t[n - 1]
        else:
            return t[k]

    H1 = (x - node(k - 1)) / (node(k) - node(k - 1))  # upward slope
    H2 = (node(k + 1) - x) / (node(k + 1) - node(k))  # downward slope
    H = np.minimum(H1, H2)
    return np.maximum(0, H)

Piecewise linear interpolation

plinterp.py

1
2
3
4
5
6
7
8
9
def plinterp(t, y):
    """
    plinterp(t,y)

    Create a piecewise linear interpolating function for data values in `y` given at nodes
    in `t`.
    """
    n = len(t) - 1
    return lambda x: np.sum(y[k] * hatfun(x, t, k) for k in range(n + 1))

Cubic spline interpolation

spinterp.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
def spinterp(t, y):
    """
    spinterp(t,y)

    Create a cubic not-a-knot spline interpolating function for data values in `y` given at nodes in `t`.
    """
    n = len(t) - 1
    h = [t[i + 1] - t[i] for i in range(n)]

    # Preliminary definitions.
    Z = np.zeros([n, n])
    I = np.eye(n)
    E = I[: n - 1, :]
    J = np.eye(n) + np.diag(-np.ones(n - 1), 1)
    H = np.diag(h)

    # Left endpoint interpolation:
    AL = np.hstack([I, Z, Z, Z])
    vL = y[:-1]

    # Right endpoint interpolation:
    AR = np.hstack([I, H, H**2, H**3])
    vR = y[1:]

    # Continuity of first derivative:
    A1 = E @ np.hstack([Z, J, 2 * H, 3 * H**2])
    v1 = np.zeros(n - 1)

    # Continuity of second derivative:
    A2 = E @ np.hstack([Z, Z, J, 3 * H])
    v2 = np.zeros(n - 1)

    # Not-a-knot conditions:
    nakL = np.hstack([np.zeros(3 * n), np.hstack([1, -1, np.zeros(n - 2)])])
    nakR = np.hstack([np.zeros(3 * n), np.hstack([np.zeros(n - 2), 1, -1])])

    # Assemble and solve the full system.
    A = np.vstack([AL, AR, A1, A2, nakL, nakR])
    v = np.hstack([vL, vR, v1, v2, 0, 0])
    z = solve(A, v)

    # Break the coefficients into separate vectors.
    rows = np.arange(n)
    a = z[rows]
    b = z[n + rows]
    c = z[2 * n + rows]
    d = z[3 * n + rows]
    S = [np.poly1d([d[k], c[k], b[k], a[k]]) for k in range(n)]

    # This function evaluates the spline when called with a value for x.
    def evaluate(x):
        f = np.zeros(x.shape)
        for k in range(n):
            # Evaluate this piece's cubic at the points inside it.
            index = (x >= t[k]) & (x <= t[k + 1])
            f[index] = S[k](x[index] - t[k])
        return f

    return evaluate

Fornberg’s algorithm for finite difference weights

fdweights.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
def fdweights(t, m):
    """
    fdweights(t,m)

    Return weights for the `m`th derivative of a function at zero using values at the
    nodes in vector `t`.
    """
    # This is a compact implementation, not an efficient one.

    def weight(t, m, r, k):
        # Recursion for one weight.
        # Input:
        #   t   nodes (vector)
        #   m   order of derivative sought
        #   r   number of nodes to use from t (<= length(t))
        #   k   index of node whose weight is found

        if (m < 0) or (m > r):  # undefined coeffs must be zero
            c = 0
        elif (m == 0) and (r == 0):  # base case of one-point interpolation
            c = 1
        else:  # generic recursion
            if k < r:
                c = t[r] * weight(t, m, r - 1, k) - m * weight(t, m - 1, r - 1, k)
                c = c / (t[r] - t[k])
            else:
                if r <= 1:
                    numer = 1.0
                else:
                    numer = np.prod(t[r-1] - t[:r-1])
                if r <= 0:
                    denom = 1.0
                else:
                    denom = np.prod(t[r] - t[:r])
                beta = numer / denom
                c = weight(t, m - 1, r - 1, r - 1) - t[r-1] * weight(t, m, r - 1, r - 1)
                c *= beta
        return c

    r = len(t) - 1
    w = np.zeros(t.shape)
    return [weight(t, m, r, k) for k in range(r + 1)]

Trapezoid formula for numerical integration

trapezoid.py

1
2
3
4
5
6
7
8
9
10
11
def trapezoid(f, a, b, n):
    """
    trapezoid(f,a,b,n)

    Apply the trapezoid integration formula for integrand `f` over interval [`a`,`b`], broken up into `n` equal pieces. Returns estimate, vector of nodes, and vector of integrand values at the nodes.
    """
    h = (b - a) / n
    t = np.linspace(a, b, n + 1)
    y = f(t)
    T = h * (np.sum(y[1:-1]) + 0.5 * (y[0] + y[-1]))
    return T, t, y

Adaptive integration

intadapt.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
def intadapt(f, a, b, tol):
    """
    intadapt(f,a,b,tol)

    Do adaptive integration to estimate the integral of `f` over [`a`,`b`] to desired
    error tolerance `tol`. Returns estimate and a vector of evaluation nodes used.
    """

    # Use error estimation and recursive bisection.
    def do_integral(a, fa, b, fb, m, fm, tol):
        # These are the two new nodes and their f-values.
        xl = (a + m) / 2
        fl = f(xl)
        xr = (m + b) / 2
        fr = f(xr)
        t = np.array([a, xl, m, xr, b])  # all 5 nodes at this level

        # Compute the trapezoid values iteratively.
        h = b - a
        T = np.zeros(3)
        T[0] = h * (fa + fb) / 2
        T[1] = T[0] / 2 + (h / 2) * fm
        T[2] = T[1] / 2 + (h / 4) * (fl + fr)

        S = (4 * T[1:] - T[:-1]) / 3  # Simpson values
        E = (S[1] - S[0]) / 15  # error estimate

        if abs(E) < tol * (1 + abs(S[1])):  # acceptable error?
            Q = S[1]  # yes--done
        else:
            # Error is too large--bisect and recurse.
            QL, tL = do_integral(a, fa, m, fm, xl, fl, tol)
            QR, tR = do_integral(m, fm, b, fb, xr, fr, tol)
            Q = QL + QR
            t = np.hstack([tL, tR[1:]])  # merge the nodes w/o duplicate
        return Q, t

    m = (b + a) / 2
    Q, t = do_integral(a, f(a), b, f(b), m, f(m), tol)
    return Q, t

Examples¶

import FNC
from numpy import *
from matplotlib.pyplot import *
from numpy.linalg import solve, norm
from scipy.interpolate import interp1d
from scipy.integrate import quad
from prettytable import PrettyTable

Section 5.1¶

Trouble in polynomial interpolation

Here are some points that we could consider to be observations of an unknown function on $[-1,1]$ .

n = 5
t = linspace(-1, 1, n + 1)
y = t**2 + t + 0.05 * sin(20 * t)
fig, ax = subplots()
plot(t, y, "o", label="data")
xlabel("$x$")
ylabel("$y$")

The polynomial interpolant, as computed using fit, looks very sensible. It’s the kind of function you’d take home to meet your parents.

p = poly1d(polyfit(t, y, n))  # interpolating polynomial
tt = linspace(-1, 1, 400)
ax.plot(tt, p(tt), label="interpolant")
ax.legend()
fig

But now consider a different set of points generated in almost exactly the same way.

n = 18
t = linspace(-1, 1, n + 1)
y = t**2 + t + 0.05 * sin(20 * t)
fig, ax = subplots()
plot(t, y, "o", label="data")
xlabel("$x$")
ylabel("$y$")

The points themselves are unremarkable. But take a look at what happens to the polynomial interpolant.

p = poly1d(polyfit(t, y, n))
ax.plot(tt, p(tt), label="interpolant")
ax.legend()
fig

Surely there must be functions that are more intuitively representative of those points!

Piecewise polynomial interpolation

Let us recall the data from Demo 5.1.1.

clf
n = 12
t = linspace(-1, 1, n + 1)
y = t**2 + t + 0.5 * sin(20 * t)
fig, ax = subplots()
scatter(t, y, label="data")
xlabel("$x$")
ylabel("$y$")

Here is an interpolant that is linear between each consecutive pair of nodes, using plinterp from Piecewise linear interpolation.

tt = linspace(-1, 1, 400)
p = interp1d(t, y, kind="linear")
ax.plot(tt, p(tt), label="piecewise linear")
ax.legend()
fig

We may prefer a smoother interpolant that is piecewise cubic:

scatter(t, y, label="data")
p = interp1d(t, y, kind="cubic")
tt = linspace(-1, 1, 400)
plot(tt, p(tt), label="cubic spline")
xlabel("$x$")
ylabel("$y$")
legend()

Conditioning of interpolation

In Demo 5.1.1 and Demo 5.1.3 we saw a big difference between polynomial interpolation and piecewise polynomial interpolation of some arbitrarily chosen data. The same effects can be seen clearly in the cardinal functions, which are closely tied to the condition numbers.

clf
n = 18
t = linspace(-1, 1, n + 1)
y = zeros(n + 1)
y[9] = 1.0
p = interp1d(t, y, kind="cubic")

scatter(t, y, label="data")
tt = linspace(-1, 1, 400)
plot(tt, p(tt), label="cardinal function")
title("Cubic spline cardinal function")
legend()

The piecewise cubic cardinal function is nowhere greater than one in absolute value. This happens to be true for all the cardinal functions, ensuring a good condition number for any interpolation with these functions. But the story for global polynomials is very different.

p = poly1d(polyfit(t, y, n))
scatter(t, y, label="data")
plot(tt, p(tt), label="cardinal function")
xlabel("$x$")
ylabel("$y$")
title("Polynomial cardinal function")
legend()

From the figure we can see that the condition number for polynomial interpolation on these nodes is at least 500.

Section 5.2¶

A look at hat functions

Let’s define a set of four nodes (i.e., $n=3$ in our formulas).

t = array([0, 0.075, 0.25, 0.55, 0.7, 1])

We plot the hat functions $H_0,\ldots,H_3$ .

x = linspace(0, 1, 300)
for k in range(6):
    plot(x, FNC.hatfun(x, t, k))
xlabel("$x$")
ylabel("$H_k(x)$")
title("Hat functions")

Using piecewise linear interpolation

We generate a piecewise linear interpolant of $f(x)=e^{\sin 7x}$ .

f = lambda x: exp(sin(7 * x))
x = linspace(0, 1, 400)
fig, ax = subplots()
plot(x, f(x), label="function")
xlabel("$x$")
ylabel("$f(x)$")

First we sample the function to create the data.

t = array([0, 0.075, 0.25, 0.55, 0.7, 1])  # nodes
y = f(t)  # function values

ax.plot(t, y, "o", label="nodes")
ax.legend()
fig

Now we create a callable function that will evaluate the piecewise linear interpolant at any $x$ , and then plot it.

p = FNC.plinterp(t, y)
ax.plot(x, p(x), label="interpolant")
ax.legend()
fig

Convergence of piecewise linear interpolation

We measure the convergence rate for piecewise linear interpolation of $e^{\sin 7x}$ over $x \in [0,1]$ .

f = lambda x: exp(sin(7 * x))
x = linspace(0, 1, 10000)  # sample the difference at many points
N = 2 ** arange(3, 11)
err = zeros(N.size)
for i, n in enumerate(N):
    t = linspace(0, 1, n + 1)  # interpolation nodes
    p = FNC.plinterp(t, f(t))
    err[i] = max(abs(f(x) - p(x)))
print(err)

As predicted, a factor of 10 in $n$ produces a factor of 100 in the error. In a convergence plot, it is traditional to have $h$ decrease from left to right, so we expect a straight line of slope -2 on a log-log plot.

order2 = 0.1 * (N / N[0]) ** (-2)
loglog(N, err, "-o", label="observed error")
loglog(N, order2, "--", label="2nd order")
xlabel("$n$")
ylabel("$\|f-p\|_\infty$")
legend()

Section 5.3¶

Cubic splines

For illustration, here is a spline interpolant using just a few nodes.

f = lambda x: exp(sin(7 * x))

x = linspace(0, 1, 500)
fig, ax = subplots()
ax.plot(x, f(x), label="function")

t = array([0, 0.075, 0.25, 0.55, 0.7, 1])  # nodes
y = f(t)  # values at nodes

xlabel("$x$")
ylabel("$y$")
ax.scatter(t, y, label="nodes")

S = FNC.spinterp(t, y)
ax.plot(x, S(x), label="spline")
ax.legend()
fig

Now we look at the convergence rate as the number of nodes increases.

N = floor(2 ** linspace(3, 8, 17)).astype(int)
err = zeros(N.size)
for i, n in enumerate(N):
    t = linspace(0, 1, n + 1)  # interpolation nodes
    p = FNC.spinterp(t, f(t))
    err[i] = max(abs(f(x) - p(x)))
print(err)

Since we expect convergence that is $O(h^4)=O(n^{-4})$ , we use a log-log graph of error and expect a straight line of slope -4.

order4 = (N / N[0]) ** (-4)
loglog(N, err, "-o", label="observed error")
loglog(N, order4, "--", label="4th order")
xlabel("$n$")
ylabel("$\|f-S\|_\infty$")
legend()

Section 5.4¶

Finite differences

If $f(x)=e^{\,\sin(x)}$ , then $f'(0)=1$ .

f = lambda x: exp(sin(x))

Here are the first two centered differences from Table 5.4.1.

h = 0.05
CD2 = (-f(-h) + f(h)) / (2*h)
CD4 = (f(-2*h) - 8*f(-h) + 8*f(h) - f(2*h)) / (12*h)
print(f"CD2 is {CD2:.9f} and CD4 is {CD4:.9f}")

Here are the first two forward differences from Table 5.4.2.

FD1 = (-f(0) + f(h)) / h
FD2 = (-3*f(0) + 4*f(h) - f(2*h)) / (2*h)
print(f"FD1 is {FD1:.9f} and FD2 is {FD2:.9f}")

Finally, here are the backward differences that come from reverse-negating the forward differences.

BD1 = (-f(-h) + f(0)) / h
BD2 = (f(-2*h) - 4*f(-h) + 3*f(0)) / (2*h)
print(f"BD1 is {BD1:.9f} and BD2 is {BD2:.9f}")

Finite differences for

f''

If $f(x)=e^{\,\sin(x)}$ , then $f''(0)=1$ .

f = lambda x: exp(sin(x))

Here is a centered estimate given by (5.4.12).

h = 0.05
CD2 = (f(-h) - 2*f(0) + f(h)) / h**2
print(f"CD2 is {CD2:.9f}")

For the same $h$ , here are forward estimates given by (5.4.13) and (5.4.14).

FD1 = (f(0) - 2*f(h) + f(2*h)) / h**2
FD2 = (2*f(0) - 5*f(h) + 4*f(2*h) - f(3*h)) / h**2
print(f"FD1 is {FD1:.9f} and FD2 is {FD2:.9f}")

Finally, here are the backward estimates that come from reversing (5.4.13) and (5.4.14).

BD1 = (f(-2*h) - 2*f(-h) + f(0)) / h**2
BD2 = (-f(-3*h) + 4*f(-2*h) - 5*f(-h) + 2*f(0)) / h**2
print(f"BD1 is {BD1:.9f} and BD2 is {BD2:.9f}")

Finite differences at arbitrary nodes

We will estimate the derivative of $\cos(x^2)$ at $x=0.5$ using five nodes.

t = array([0.35, 0.5, 0.57, 0.6, 0.75])   # nodes
f = lambda x: cos(x**2)
dfdx = lambda x: -2 * x * sin(x**2)
exact_value = dfdx(0.5)

We have to shift the nodes so that the point of estimation for the derivative is at $x=0$ . (To subtract a scalar from a vector, we must use the .- operator.)

w = FNC.fdweights(t - 0.5, 1)

The finite-difference formula is a dot product (i.e., inner product) between the vector of weights and the vector of function values at the nodes.

fd_value = dot(w, f(t))

We can reproduce the weights in the finite-difference tables by using equally spaced nodes with $h=1$ . For example, here is a one-sided formula at four nodes.

print(FNC.fdweights(linspace(0, 3, 4), 1))

Section 5.5¶

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)

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)

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()

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)

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).

Section 5.6¶

Numerical integration

The antiderivative of $e^x$ is, of course, itself. That makes evaluation of $\int_0^1 e^x\,dx$ by the Fundamental Theorem trivial.

exact = exp(1) - 1

The module scipy.integrate has multiple functions that estimate the value of an integral numerically without finding the antiderivative first. As you can see here, it’s often just as accurate.

Q, errest = quad(exp, 0, 1, epsabs=1e-13, epsrel=1e-13)
print(Q)

The numerical approach is also far more robust. For example, $e^{\,\sin x}$ has no useful antiderivative. But numerically, it’s no more difficult.

Q, errest = quad(lambda x: exp(sin(x)), 0, 1, epsabs=1e-13, epsrel=1e-13)
print(Q)

When you look at the graphs of these functions, what’s remarkable is that one of these areas is basic calculus while the other is almost impenetrable analytically. From a numerical standpoint, they are practically the same problem.

x = linspace(0, 1, 300)
subplot(1, 2, 1)
plot(x, exp(x))
ylim([0, 2.7]), title("exp(x)")
subplot(1, 2, 2)
plot(x, exp(sin(x)))
ylim([0, 2.7]), title("exp(sin(x))");

Trapezoid integration

We will approximate the integral of the function $f(x)=e^{\sin 7x}$ over the interval $[0,2]$ .

f = lambda x: exp(sin(7 * x))
a, b = 0, 2

In lieu of the exact value, we will use the quad function to find an accurate result.

I, errest = quad(f, a, b, epsabs=1e-13, epsrel=1e-13)
print(f"Integral = {I:.14f}")

Here is the trapezoid result at $n=40$ , and its error.

T, t, y = FNC.trapezoid(f, a, b, 40)
print(f"Trapezoid estimate is {T:.14f} with error {I - T:.2e}")

In order to check the order of accuracy, we increase $n$ by orders of magnitude and observe how the error decreases.

n_ = 40 * 2 ** arange(6)
err = zeros(size(n_))
print("     n     error")
for k, n in enumerate(n_):
    T, t, y = FNC.trapezoid(f, a, b, n)
    err[k] = I - T
    print(f"{n:6d}   {err[k]:8.3e} ")

Each increase by a factor of 10 in $n$ cuts the error by a factor of about 100, which is consistent with second-order convergence. Another check is that a log-log graph should give a line of slope -2 as $n\to\infty$ .

loglog(n_, abs(err), "-o", label="results")
loglog(n_, 3e-3 * (n_ / n_[0]) ** (-2), "--", label="2nd order")
gca().invert_xaxis()
xlabel("$n$")
ylabel("error")
legend()
title("Convergence of trapezoidal integration");

Integration by extrapolation

We estimate $\displaystyle\int_0^2 x^2 e^{-2x}\, dx$ using extrapolation. First we use quadgk to get an accurate value.

f = lambda x: x**2 * exp(-2 * x)
a = 0
b = 2
I, errest = quad(f, a, b, epsabs=1e-13, epsrel=1e-13)
print(f"Integral = {I:.14f}")

We start with the trapezoid formula on $n=N$ nodes.

N = 20    # the coarsest formula
n = N
h = (b - a) / n
t = h * arange(n + 1)
y = f(t)

We can now apply weights to get the estimate $T_f(N)$ .

T = zeros(3)
T[0] = h * (sum(y[1:-1]) + y[0] / 2 + y[-1] / 2)
print(f"error (2nd order): {I - T[0]:.2e}")

Now we double to $n=2N$ , but we only need to evaluate $f$ at every other interior node and apply (5.6.18).

n = 2 * n
h = h / 2
t = h * arange(n + 1)
T[1] = T[0] / 2 + h * sum(f(t[1:-1:2]))
print("error (2nd order):", I - T[:2])

As expected for a second-order estimate, the error went down by a factor of about 4. We can repeat the same code to double $n$ again.

n = 2 * n
h = h / 2
t = h * arange(n + 1)
T[2] = T[1] / 2 + h * sum(f(t[1:-1:2]))
print("error (2nd order):", I - T[:3])

Let us now do the first level of extrapolation to get results from Simpson’s formula. We combine the elements T[i] and T[i+1] the same way for $i=1$ and $i=2$ .

S = array([(4 * T[i + 1] - T[i]) / 3 for i in range(2)])
print("error (4th order):", I - S)

With the two Simpson values $S_f(N)$ and $S_f(2N)$ in hand, we can do one more level of extrapolation to get a sixth-order accurate result.

R = (16 * S[1] - S[0]) / 15
print("error (6th order):", I - R)

We can make a triangular table of the errors:

err = nan * ones((3, 3))
err[0, :] = I - T
err[1, 1:] = I - S
err[2, 2] = I - R
results = PrettyTable(["2nd order", "4th order", "6th order"])
results.add_rows(err.T)
print(results)

If we consider the computational time to be dominated by evaluations of $f$ , then we have obtained a result with about twice as many accurate digits as the best trapezoid result, at virtually no extra cost.

Section 5.7¶

Motivation for adaptive integration

This function gets increasingly oscillatory as $x$ increases.

f = lambda x: (x + 1) ** 2 * cos((2 * x + 1) / (x - 4.3))
x = linspace(0, 4, 600)
plot(x, f(x))
xlabel("$x$")
ylabel("$f(x)$");

Accordingly, the trapezoid rule is more accurate on the left half of this interval than on the right half.

n_ = 50 * 2 ** arange(4)
Tleft = zeros(4)
Tright = zeros(4)
for i, n in enumerate(n_):
    Tleft[i] = FNC.trapezoid(f, 0, 2, n)[0]
    Tright[i] = FNC.trapezoid(f, 2, 4, n)[0]
print("left half:", Tleft)
print("right half:", Tright)

left_val, err = quad(f, 0, 2, epsabs=1e-13, epsrel=1e-13)
right_val, err = quad(f, 2, 4, epsabs=1e-13, epsrel=1e-13)

print("    n     left error   right error")
for k in range(n_.size):
    print(f"  {n_[k]:4}    {Tleft[k]-left_val:8.3e}    {Tright[k]-right_val:8.3e}")

Both the picture and the numerical results suggest that more nodes should be used on the right half of the interval than on the left half.

Using adaptive integration

We’ll integrate the function from Demo 5.7.1.

f = lambda x: (x + 1) ** 2 * cos((2 * x + 1) / (x - 4.3))
I, errest = quad(f, 0, 4, epsabs=1e-12, epsrel=1e-12)
print("integral:", I)    # 'exact' value

We perform the integration and show the nodes selected underneath the curve.

Q, t = FNC.intadapt(f, 0, 4, 0.001)
print("number of nodes:", t.size)

x = linspace(0, 4, 600)
plot(x, f(x), "k")
stem(t, f(t))
xlabel("$x$"); ylabel("$f(x)$");

The error turns out to be a bit more than we requested. It’s only an estimate, not a guarantee.

print("error:", I - Q)

Let’s see how the number of integrand evaluations and the error vary with the requested tolerance.

tol_ = 10.0 ** arange(-4, -12, -1)
err_ = zeros(tol_.size)
num_ = zeros(tol_.size, dtype=int)
print("    tol         error     # f-evals")
for i, tol in enumerate(tol_):
    Q, t = FNC.intadapt(f, 0, 4, tol)
    err_[i] = I - Q
    num_[i] = t.size
    print(f"  {tol:6.1e}    {err_[i]:10.3e}    {num_[i]:6d}")

As you can see, even though the errors are not smaller than the estimates, the two columns decrease in tandem. If we consider now the convergence not in $h$ , which is poorly defined now, but in the number of nodes actually chosen, we come close to the fourth-order accuracy of the underlying Simpson scheme.

loglog(num_, abs(err_), "-o", label="results")
order4 = 0.01 * (num_ / num_[0]) ** (-4)
loglog(num_, order4, "--", label="$O(n^{-4})$")
xlabel("number of nodes"), ylabel("error")
legend()
title("Convergence of adaptive quadrature");