Adaptive integration Tobin A. Driscoll Richard J. Braun
To this point, we have used only equally spaced nodes to compute integrals. Yet there are problems in which non-uniformly distributed nodes would clearly be more appropriate, as demonstrated in Example 5.7.1
Example 5.7.1 (Motivation for adaptive integration)
Example 5.7.1 This function gets increasingly oscillatory as x x x
using Plots
f = x -> (x + 1)^2 * cos((2x + 1) / (x - 4.3))
plot(f, 0, 4, xlabel=L"x", ylabel=L"f(x)")Accordingly, the trapezoid rule is more accurate on the left half of this interval than on the right half.
using QuadGK
left_val, _ = quadgk(f, 0, 2, atol=1e-14, rtol=1e-14)
right_val, _ = quadgk(f, 2, 4, atol=1e-14, rtol=1e-14)
n = [50 * 2^k for k in 0:3]
err = zeros(length(n), 2)
for (k, n) in enumerate(n)
T, _ = FNC.trapezoid(f, 0, 2, n)
err[k, 1] = T - left_val
T, _ = FNC.trapezoid(f, 2, 4, n)
err[k, 2] = T - right_val
end
@pt :header=["n", "left error", "right error"] [n err]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.
Example 5.7.1 This function gets increasingly oscillatory as x x x
f = @(x) (x + 1).^2 .* cos((2 * x + 1) ./ (x - 4.3));
clf
fplot(f, [0, 4], 2000)
xlabel('x'), ylabel(('f(x)'));Accordingly, the trapezoid rule is more accurate on the left half of this interval than on the right half.
left_val = integral(f, 0, 2, abstol=1e-14, reltol=1e-14);
right_val = integral(f, 2, 4, abstol=1e-14, reltol=1e-14);
n = round(50 * 2 .^ (0:3)');
err = zeros(length(n), 2);
for i = 1:length(n)
T = trapezoid(f, 0, 2, n(i));
err(i, 1) = T - left_val;
T = trapezoid(f, 2, 4, n(i));
err(i, 2) = T - right_val;
end
table(n, err(:, 1), err(:, 2), variableNames=["n", "left error", "right error"])Unrecognized function or variable 'trapezoid'. 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.
Example 5.7.1 This function gets increasingly oscillatory as x x x
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 half: [2.00419122 2.00605957 2.00652661 2.00664337]
right half: [-4.32798637 -4.73621129 -4.80966839 -4.82666144]
from scipy.integrate import quad
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}") n left error right error
50 -2.491e-03 5.042e-01
100 -6.227e-04 9.600e-02
200 -1.557e-04 2.255e-02
400 -3.892e-05 5.554e-03
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.
We would like an algorithm that automatically detects and reacts to a situation like that in Example 5.7.1 adaptivity .
5.7.1 Error estimation ¶ Ideally, we would like to make adaptation decisions based on the error of the integration result. Knowing the error exactly would be equivalent to knowing the exact answer, but we can estimate it using the extrapolation technique of Numerical integration . Consider the Simpson formula (5.6.15)
We expect this method to be fourth-order accurate, i.e.,
∫ a b f ( x ) d x = S f ( 2 n ) + O ( n − 4 ) , \int_a^b f(x)\, dx = S_f(2n) + O(n^{-4}), ∫ a b f ( x ) d x = S f ( 2 n ) + O ( n − 4 ) , We can further extrapolate to sixth-order accuracy using (5.6.17)
By virtue of higher order of accuracy, R f ( 4 n ) R_f(4n) R f ( 4 n ) S f ( 4 n ) S_f(4n) S f ( 4 n )
E = R f ( 4 n ) − S f ( 4 n ) = S f ( 4 n ) − S f ( 2 n ) 15 . E = R_f(4n) - S_f(4n) = \frac{S_f(4n) - S_f(2n)}{15}. E = R f ( 4 n ) − S f ( 4 n ) = 15 S f ( 4 n ) − S f ( 2 n ) . 5.7.2 Divide and conquer ¶ If ∣ E ∣ |E| ∣ E ∣ 1020 . Requiring ∣ E ∣ < δ ≪ 1 |E| < \delta\ll 1 ∣ E ∣ < δ ≪ 1
∫ 1 0 − 6 2 π 2 sin x d x ≈ − 1 0 − 12 . \int_{10^{-6}}^{2\pi} 2 \sin x\, dx \approx -10^{-12}. ∫ 1 0 − 6 2 π 2 sin x d x ≈ − 1 0 − 12 . We are likely to sample values of the integrand that are larger than, say, 1 / 2 1/2 1/2 O ( 1 ) O(1) O ( 1 )
Typically, we use both relative and absolute error, stopping when either one is considered small enough. Algebraically, the test is
∣ E ∣ < δ a + δ r ∣ S f ( n ) ∣ , |E| < \delta_a + \delta_r |S_f(n)|, ∣ E ∣ < δ a + δ r ∣ S f ( n ) ∣ , where δ a \delta_a δ a δ r \delta_r δ r
When ∣ E ∣ |E| ∣ E ∣ (5.7.6) [ a , b ] [a,b] [ a , b ]
∫ a b f ( x ) d x = ∫ a ( a + b ) / 2 f ( x ) d x + ∫ ( a + b ) / 2 b f ( x ) d x , \int_a^b f(x)\, dx = \int_a^{(a+b)/2} f(x)\, dx + \int_{(a+b)/2}^b f(x)\, dx, ∫ a b f ( x ) d x = ∫ a ( a + b ) /2 f ( x ) d x + ∫ ( a + b ) /2 b f ( x ) d x , and independently compute estimates to each of the half-length integrals. Each of these half-sized computations recursively applies Simpson’s formula and the error estimation criterion, making further bisections as necessary. Such an approach is called divide and conquer in computer science: recursively split the problem into easier pieces and glue the results together.
5.7.3 Implementation ¶ It is typical to use just the minimal formula S f ( 4 ) S_f(4) S f ( 4 ) E E E S f ( 4 ) S_f(4) S f ( 4 ) T f ( 1 ) T_f(1) T f ( 1 ) T f ( 2 ) T_f(2) T f ( 2 ) T f ( 4 ) T_f(4) T f ( 4 ) (5.6.18) Example 5.6.3 T f ( 4 ) T_f(4) T f ( 4 )
Function 5.7.1 n = 1 n=1 n = 1 n = 2 n=2 n = 2 n = 4 n=4 n = 4 (5.6.18) (5.7.4)
If the error estimate passes the test (5.7.6)
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
"""
intadapt(f, a, b, tol)
Adaptively integrate `f` over [`a`,`b`] to within target error
tolerance `tol`. Returns the estimate and a vector of evaluation
nodes.
"""
function intadapt(f, a, b, tol, fa = f(a), fb = f(b), m = (a + b) / 2, fm = f(m))
# Use error estimation and recursive bisection.
# These are the two new nodes and their f-values.
xl = (a + m) / 2
fl = f(xl)
xr = (m + b) / 2
fr = f(xr)
# Compute the trapezoid values iteratively.
h = (b - a)
T = [0.0, 0.0, 0.0]
T[1] = h * (fa + fb) / 2
T[2] = T[1] / 2 + (h / 2) * fm
T[3] = T[2] / 2 + (h / 4) * (fl + fr)
S = (4T[2:3] - T[1:2]) / 3 # Simpson values
E = (S[2] - S[1]) / 15 # error estimate
if abs(E) < tol * (1 + abs(S[2])) # acceptable error?
Q = S[2] # yes--done
nodes = [a, xl, m, xr, b] # all nodes at this level
else
# Error is too large--bisect and recurse.
QL, tL = intadapt(f, a, m, tol, fa, fm, xl, fl)
QR, tR = intadapt(f, m, b, tol, fm, fb, xr, fr)
Q = QL + QR
nodes = [tL; tR[2:end]] # merge the nodes w/o duplicate
end
return Q, nodes
endAbout the code
The intended way for a user to call Function 5.7.1 f, a, b, and tol provided. We then use default values on the other parameters to compute the function values at the endpoints, the interval’s midpoint, and the function value at the midpoint. Recursive calls from within the function itself will provide all of that information, since it was already calculated along the way.
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
function [Q,t] = intadapt(f,a,b,tol)
%INTADAPT Adaptive integration with error estimation.
% Input:
% f integrand (function)
% a,b interval of integration (scalars)
% tol acceptable error
% Output:
% Q approximation to integral(f, a, b)
% t vector of nodes used
m = (b + a) / 2;
[Q, t] = do_integral(a, f(a), b, f(b), m, f(m), tol);
% Use error estimation and recursive bisection.
function [Q, t] = 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 = [a; xl; m; xr; b]; % all 5 nodes at this level
% Compute the trapezoid values iteratively.
h = b - a;
T(1) = h * (fa + fb) / 2;
T(2) = T(1)/2 + (h / 2) * fm;
T(3) = T(2)/2 + (h / 4) * (fl + fr);
S = (4 * T(2:3) - T(1:2)) / 3; % Simpson values
E = (S(2) - S(1)) / 15; % error estimate
if abs(E) < tol * (1 + abs(S(2))) % acceptable error?
Q = S(2); % 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 = [tL; tR(2:end)]; % merge the nodes w/o duplicate
end
end
end % main functionAbout the code
The intended way for a user to call Function 5.7.1 f, a, b, and tol provided. We then use default values on the other parameters to compute the function values at the endpoints, the interval’s midpoint, and the function value at the midpoint. Recursive calls from within the function itself will provide all of that information, since it was already calculated along the way.
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
About the code
The intended way for a user to call Function 5.7.1 f, a, b, and tol provided. We then use default values on the other parameters to compute the function values at the endpoints, the interval’s midpoint, and the function value at the midpoint. Recursive calls from within the function itself will provide all of that information, since it was already calculated along the way.
Example 5.7.2 (Using adaptive integration)
Example 5.7.2 We’ll integrate the function from Example 5.7.1
f = x -> (x + 1)^2 * cos((2x + 1) / (x - 4.3));
We perform the integration and show the nodes selected underneath the curve.
A, t = FNC.intadapt(f, 0, 4, 0.001)
@show num_nodes = length(t);
plot(f, 0, 4;
color=:black, legend=:none,
xlabel=L"x", ylabel=L"f(x)",
title="Adaptive node selection")
plot!(t, f.(t), seriestype=:sticks, m=(:o, 2))The error turns out to be a bit more than we requested. It’s only an estimate, not a guarantee.
Q, _ = quadgk(f, 0, 4, atol=1e-14, rtol=1e-14); # 'exact' value
println("error: $(Q-A)");error: -0.022002813037627522
Let’s see how the number of integrand evaluations and the error vary with the requested tolerance.
tol = [1 / 10^k for k in 4:14]
err, n = [], []
for tol in 10.0 .^ (-4:-1:-14)
A, t = FNC.intadapt(f, 0, 4, tol)
push!(err, Q - A)
push!(n, length(t))
end
@pt :header=["tolerance", "error", "number of nodes"] [tol err n][1:2:end, :]As you can see, even though the errors are not smaller than the tolerances, the two columns decrease in tandem. If we consider now the convergence not in h h h
plot(n, abs.(err);
m=:o, label="results",
xaxis=(:log10, "number of nodes"), yaxis=(:log10, "error"),
title="Convergence of adaptive integration")
order4 = @. 0.01 * (n / n[1])^(-4)
plot!(n, order4, l=:dash, label=L"O(n^{-4})")Example 5.7.2 We’ll integrate the function from Example 5.7.1
f = @(x) (x + 1).^2 .* cos((2 * x + 1) ./ (x - 4.3));
We perform the integration and show the nodes selected underneath the curve.
[Q, t] = intadapt(f, 0, 4, 0.001);
clf, fplot(f, [0, 4], 2000)
hold on
stem(t, f(t), '.-')
title('Adaptive node selection')
xlabel('x'), ylabel('f(x)')
fprintf("number of nodes = %d", length(t))The error turns out to be a bit more than we requested. It’s only an estimate, not a guarantee.
I = integral(f, 0, 4, abstol=1e-14, reltol=1e-14); % 'exact' value
fprintf("error = %.2e", abs(Q - I))Let’s see how the number of integrand evaluations and the error vary with the requested tolerance.
tol = 1 ./ 10.^(4:14)';
err = zeros(size(tol));
n = zeros(size(tol));
for k = 1:length(tol)
[A, t] = intadapt(f, 0, 4, tol(k));
err(k) = I - A;
n(k) = length(t);
end
table(tol, err, n, variableNames=["tolerance", "error", "number of nodes"])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 h h
clf
loglog(n, abs(err), "-o", displayname="results")
xlabel("number of nodes"), ylabel("error")
title("Convergence of adaptive integration")
order4 = 0.1 * abs(err(end)) * (n / n(end)).^(-4);
hold on
loglog(n, order4, "k--", displayname="O(n^{-4})")
legend();Example 5.7.2 We’ll integrate the function from Example 5.7.1
from scipy.integrate import quad
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' valueintegral: -2.8255333734374504
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.
error: -0.02200281303763152
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}") tol error # f-evals
1.0e-04 -4.195e-04 113
1.0e-05 4.790e-05 181
1.0e-06 6.314e-06 297
1.0e-07 -6.639e-07 489
1.0e-08 7.181e-08 757
1.0e-09 1.265e-08 1193
1.0e-10 -8.441e-10 2009
1.0e-11 2.612e-11 3157
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 h h
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");Although adaptivity and the error estimation that goes with it can be very powerful, they come at some cost. The error estimation cannot be universally perfect, so sometimes the answer will not be as accurate as requested, and sometimes the function will be evaluated more times than necessary. Subtle problems may arise when the integral is a step within a larger computation (see Exercise 5.7.6
5.7.4 Exercises ¶ ⌨ For each integral below, use Function 5.7.1 1 0 − 2 , 1 0 − 3 , … , 1 0 − 12 10^{-2},10^{-3},\ldots,10^{-12} 1 0 − 2 , 1 0 − 3 , … , 1 0 − 12 Example 5.7.2 Bailey et al. (2005)
(a) ∫ 0 1 x log ( 1 + x ) d x = 1 4 \displaystyle \int_0^1 x\log(1+x)\, dx = \frac{1}{4} ∫ 0 1 x log ( 1 + x ) d x = 4 1
(b) ∫ 0 1 x 2 tan − 1 x d x = π − 2 + 2 log 2 12 \displaystyle \int_0^1 x^2 \tan^{-1}x\, dx = \frac{\pi-2+2\log 2}{12} ∫ 0 1 x 2 tan − 1 x d x = 12 π − 2 + 2 log 2
(c) ∫ 0 π / 2 e x cos x d x = e π / 2 − 1 2 \displaystyle \int_0^{\pi/2}e^x \cos x\, dx = \frac{e^{\pi/2}-1}{2} ∫ 0 π /2 e x cos x d x = 2 e π /2 − 1
(d) ∫ 0 1 x log ( x ) d x = − 4 9 \displaystyle \int_{0}^1 \sqrt{x} \log(x) \, dx = -\frac{4}{9} ∫ 0 1 x log ( x ) d x = − 9 4 x → 0 x\to 0 x → 0 f ( 0 ) f(0) f ( 0 ) x = ϵ mach x=\macheps x = ϵ mach
(e) ∫ 0 1 1 − x 2 d x = π 4 \displaystyle \int_0^1 \sqrt{1-x^2}\, dx = \frac{\pi}{4} ∫ 0 1 1 − x 2 d x = 4 π
⌨ For each integral below: (i) use quadgk to find the value to at least 12 digits; (ii) use Function 5.7.1 10-8 ; (iii) compute the absolute error and the number of nodes used; (iv) use the O ( h 2 ) O(h^2) O ( h 2 ) (5.6.9) 10-8 .
(a) ∫ 0.1 3 sech ( sin ( 1 / x ) ) d x \displaystyle \int_{0.1}^3 \operatorname{sech}(\sin(1/x))\, d x ∫ 0.1 3 sech ( sin ( 1/ x )) d x
(b) ∫ − 0.9 9 ln ( ( x + 1 ) 3 ) ) d x \rule[2em]{0pt}{0pt} \displaystyle\int_{-0.9}^9 \ln((x+1)^3))\, d x ∫ − 0.9 9 ln (( x + 1 ) 3 )) d x
(c) ∫ − π π cos ( x 3 ) d x \rule[2em]{0pt}{0pt} \displaystyle\int_{-\pi}^\pi \cos(x^3)\, d x ∫ − π π cos ( x 3 ) d x
⌨ An integral such as ∫ 0 1 x − γ d x \displaystyle \int_0^1 x^{-\gamma}\, dx ∫ 0 1 x − γ d x γ > 0 \gamma>0 γ > 0 improper integral. It has a finite value if γ < 1 \gamma<1 γ < 1 f ( t 0 ) f(t_0) f ( t 0 )
Using Function 5.7.1 γ = 2 / 3 \gamma=2/3 γ = 2/3 ϵ = 1 0 − 15 , 1 0 − 16 , … , 1 0 − 45 \epsilon=10^{-15},10^{-16},\ldots,10^{-45} ϵ = 1 0 − 15 , 1 0 − 16 , … , 1 0 − 45
⌨ A curious consequence of our logic in Function 5.7.1 S f ( 2 n ) S_f(2n) S f ( 2 n ) R f ( 4 n ) R_f(4n) R f ( 4 n ) local extrapolation .
Modify Function 5.7.1 Exercise 5.7.1
⌨ The sine integral function is defined by
Si ( x ) = ∫ 0 x sin z z d z . \operatorname{Si}(x) = \int_0^x \frac{\sin z}{z}\, dz. Si ( x ) = ∫ 0 x z sin z d z . Use Function 5.7.1 [ 1 , 10 ] [1,10] [ 1 , 10 ] ϵ mach \macheps ϵ mach
⌨ Adaptive integration can have subtle drawbacks. This exercise is based on the error function , a smooth function defined as
erf ( x ) = 2 π ∫ 0 x e − s 2 d s . \operatorname{erf}(x) = \frac{2}{\sqrt{\pi}}\int_0^x e^{-s^2}\,ds. erf ( x ) = π 2 ∫ 0 x e − s 2 d s . (a) Define a function g g g Function 5.6.1 n = 300 n=300 n = 300 g ( x ) − erf ( x ) g(x)-\operatorname{erf}(x) g ( x ) − erf ( x ) [ 0 , 3 ] [0,3] [ 0 , 3 ]
(b) Define another approximation h h h Function 5.7.1 10-7 . Plot the error in h h h
(c) Suppose you wished to find x x x erf ( x ) = . 95 \operatorname{erf}(x) = .95 erf ( x ) = .95
Bailey, D. H., Jeyabalan, K., & Li, X. S. (2005). A Comparison of Three High-Precision Quadrature Schemes. Experimental Mathematics , 14 (3), 317–329. 10.1080/10586458.2005.10128931