Lecture Materials
Adaptive methods#
Error control#
It is important to control the error of numerical integration since we are often interested in the precision of the result.
Let us consider the error estimate of the rectangle/trapezoidal method. For rectangle/trapezoidal method we know that the error scales with \(h\) as \(\varepsilon = c h^2\). Let us double the number of steps. We have \(h_2 = h_1 / 2\). Then, \(\varepsilon_2 = I - I_2 = c h_2^2\), and \(\varepsilon_1 = I - I_1 = 4 c h_2^2\). Therefore, \(\varepsilon_2 \simeq (I_2 - I_1) / 3\).
More generally, the error estimate of \(k\)th iteration is
\[
\varepsilon_k \simeq (I_k - I_{k-1}) / 3.
\]
We can thus continue to double the number of subintervals until the desired precision is reached.
Let us implement this in Python.
Adaptive rectangle rule#
# Rectangle rule for numerical integration with adaptive step
def rectangle_rule_adaptive(f, a, b, nst = 1, tol = 1.e-8, max_iterations = 16):
Iprev = 0.
n = nst
Iprev = rectangle_rule(f, a, b, n)
print("Iteration: {0:5}, I = {1:20.15f}".format(1, Iprev))
for k in range(1, max_iterations):
n *= 2
Inew = rectangle_rule(f, a, b, n)
ek = (Inew - Iprev) / 3.
print("Iteration: {0:5}, I = {1:20.15f}, error estimate = {2:10.15f}".format(k+1, Inew, ek))
if (abs(ek) < tol):
return Inew
Iprev = Inew
print("Failed to achieve the desired accuracy after", max_iterations,"iterations")
return Inew
Let illustrate this with an example.
# The function example we will use
# Overwrite as applicable
flabel = 'x^4 - 2x + 2'
def f(x):
return x**4 - 2*x + 2
flimit_a = 0.
flimit_b = 2.
a = flimit_a
b = flimit_b
print("Computing the integral of",flabel, "over the interval (",a,",",b,") using adaptive rectangle rule")
res = rectangle_rule_adaptive(f,a,b)
Computing the integral of x^4 - 2x + 2 over the interval ( 0.0 , 2.0 ) using adaptive rectangle rule
Iteration: 1, I = 2.000000000000000
Iteration: 2, I = 5.125000000000000, error estimate = 1.041666666666667
Iteration: 3, I = 6.070312500000000, error estimate = 0.315104166666667
Iteration: 4, I = 6.316894531250000, error estimate = 0.082194010416667
Iteration: 5, I = 6.379180908203125, error estimate = 0.020762125651042
Iteration: 6, I = 6.394792556762695, error estimate = 0.005203882853190
Iteration: 7, I = 6.398697972297668, error estimate = 0.001301805178324
Iteration: 8, I = 6.399674482643604, error estimate = 0.000325503448645
Iteration: 9, I = 6.399918620008975, error estimate = 0.000081379121790
Iteration: 10, I = 6.399979654961498, error estimate = 0.000020344984174
Iteration: 11, I = 6.399994913737828, error estimate = 0.000005086258777
Iteration: 12, I = 6.399998728434201, error estimate = 0.000001271565458
Iteration: 13, I = 6.399999682108611, error estimate = 0.000000317891470
Iteration: 14, I = 6.399999920527143, error estimate = 0.000000079472844
Iteration: 15, I = 6.399999980131772, error estimate = 0.000000019868210
Iteration: 16, I = 6.399999995032923, error estimate = 0.000000004967050
Adaptive trapezoidal rule#
# Trapezoidal rule for numerical integration with adaptive step
def trapezoidal_rule_adaptive(f, a, b, nst = 1, tol = 1.e-8, max_iterations = 16):
Iprev = 0.
n = nst
Iprev = trapezoidal_rule(f, a, b, n)
print("Iteration: {0:5}, I = {1:20.15f}".format(1, Iprev))
for k in range(1, max_iterations):
n *= 2
Inew = trapezoidal_rule(f, a, b, n)
ek = (Inew - Iprev) / 3.
print("Iteration: {0:5}, I = {1:20.15f}, error estimate = {2:10.15f}".format(k+1, Inew, ek))
if (abs(ek) < tol):
return Inew
Iprev = Inew
print("Failed to achieve the desired accuracy after", max_iterations,"iterations")
return Inew
a = flimit_a
b = flimit_b
print("Computing the integral of",flabel, "over the interval (",a,",",b,") using adaptive trapezoidal rule")
res =trapezoidal_rule_adaptive(f,a,b)
Computing the integral of x^4 - 2x + 2 over the interval ( 0.0 , 2.0 ) using adaptive trapezoidal rule
Iteration: 1, I = 16.000000000000000
Iteration: 2, I = 9.000000000000000, error estimate = -2.333333333333333
Iteration: 3, I = 7.062500000000000, error estimate = -0.645833333333333
Iteration: 4, I = 6.566406250000000, error estimate = -0.165364583333333
Iteration: 5, I = 6.441650390625000, error estimate = -0.041585286458333
Iteration: 6, I = 6.410415649414062, error estimate = -0.010411580403646
Iteration: 7, I = 6.402604103088379, error estimate = -0.002603848775228
Iteration: 8, I = 6.400651037693024, error estimate = -0.000651021798452
Iteration: 9, I = 6.400162760168314, error estimate = -0.000162759174903
Iteration: 10, I = 6.400040690088645, error estimate = -0.000040690026556
Iteration: 11, I = 6.400010172525072, error estimate = -0.000010172521191
Iteration: 12, I = 6.400002543131352, error estimate = -0.000002543131240
Iteration: 13, I = 6.400000635782950, error estimate = -0.000000635782801
Iteration: 14, I = 6.400000158945742, error estimate = -0.000000158945736
Iteration: 15, I = 6.400000039736406, error estimate = -0.000000039736446
Iteration: 16, I = 6.400000009934106, error estimate = -0.000000009934100
Adaptive Simpson’s rule#
The error estimate of Simpson’s rule is \(\varepsilon = c h^4\). For this reason the error estimate of \(k\)th iteration is
\[
\varepsilon_k \simeq (I_k - I_{k-1}) / 15.
\]
# Simpson's rule for numerical integration with adaptive step
def simpson_rule_adaptive(f, a, b, nst = 2, tol = 1.e-8, max_iterations = 16):
Iprev = 0.
n = nst
Iprev = simpson_rule(f, a, b, n)
print("Iteration: {0:5}, I = {1:20.15f}".format(1, Iprev))
for k in range(1, max_iterations):
n *= 2
Inew = simpson_rule(f, a, b, n)
ek = (Inew - Iprev) / 15.
print("Iteration: {0:5}, I = {1:20.15f}, error estimate = {2:10.15f}".format(k+1, Inew, ek))
if (abs(ek) < tol):
return Inew
Iprev = Inew
print("Failed to achieve the desired accuracy after", max_iterations,"iterations")
return Inew
a = flimit_a
b = flimit_b
print("Computing the integral of",flabel, "over the interval (",a,",",b,") using adaptive Simpson's rule")
res = simpson_rule_adaptive(f,a,b,2)
Computing the integral of x^4 - 2x + 2 over the interval ( 0.0 , 2.0 ) using adaptive Simpson's rule
Iteration: 1, I = 6.666666666666666
Iteration: 2, I = 6.416666666666666, error estimate = -0.016666666666667
Iteration: 3, I = 6.401041666666666, error estimate = -0.001041666666667
Iteration: 4, I = 6.400065104166666, error estimate = -0.000065104166667
Iteration: 5, I = 6.400004069010416, error estimate = -0.000004069010417
Iteration: 6, I = 6.400000254313150, error estimate = -0.000000254313151
Iteration: 7, I = 6.400000015894571, error estimate = -0.000000015894572
Iteration: 8, I = 6.400000000993410, error estimate = -0.000000000993411
Example: Runge function#
Recall the problematic Runge function
\[
f(x) = \frac{1}{1+25x^2}
\]
Let us compute the integral of this function over the interval \([-2,2]\)
rungelabel = "Runge function"
def runge(x):
return 1./(25*x**2 + 1.)
a = -2.
b = 2.
print("Computing the integral of",rungelabel, "over the interval (",a,",",b,") using adaptive trapezoidal rule")
res = trapezoidal_rule_adaptive(runge,a,b,4,1.e-10)
Computing the integral of Runge function over the interval ( -2.0 , 2.0 ) using adaptive trapezoidal rule
Iteration: 1, I = 1.086824067022087
Iteration: 2, I = 0.698810316902099, error estimate = -0.129337916706663
Iteration: 3, I = 0.596649043819530, error estimate = -0.034053757694190
Iteration: 4, I = 0.588479663841841, error estimate = -0.002723126659230
Iteration: 5, I = 0.588444691123849, error estimate = -0.000011657572664
Iteration: 6, I = 0.588449474263155, error estimate = 0.000001594379768
Iteration: 7, I = 0.588450670842736, error estimate = 0.000000398859860
Iteration: 8, I = 0.588450970000918, error estimate = 0.000000099719394
Iteration: 9, I = 0.588451044791294, error estimate = 0.000000024930125
Iteration: 10, I = 0.588451063488940, error estimate = 0.000000006232549
Iteration: 11, I = 0.588451068163354, error estimate = 0.000000001558138
Iteration: 12, I = 0.588451069331961, error estimate = 0.000000000389535
Iteration: 13, I = 0.588451069624111, error estimate = 0.000000000097383
Romberg method#
Romberg method is a generalization of the above logic for cancelling higher-order errors. The corresponding procedure is called Richardson extrapolation and forms a quintessential method for computing integrals using equidistant grids.
Derivation#
Recall how we estimated the \(\mathcal{O}(h^2)\) error in the \(k\)th step of the trapezoidal method as
\[
I - I_k = \frac{I_k - I_{k-1}}{3} + \mathcal{O}(h^4).
\]
The integral can therefore be estimated to \(\mathcal{O}(h^4)\) order as
\[
I = I_k + \frac{I_k - I_{k-1}}{3} + \mathcal{O}(h^4),
\]
which is in fact nothing else but the Simpson rule.
We can denote \(R_{k,0} = I_k\) and \(R_{k,1} = R_{k,0} + \frac{R_{k,0} - R_{k-1,0}}{3}\). As seen above
\[
I = R_{k,1} + \mathcal{O}(h^4).
\]
Repeating this process to eliminate the \(\mathcal{O}(h^4)\) we get a higher-order estimate
\[
R_{k,2} = R_{k,1} + \frac{R_{k,1} - R_{k-1,1}}{15}
\]
which accurate to order \(\mathcal{O}(h^6)\).
The general formula for an integral estimate of order \(\mathcal{O}(h^{2m+2})\) reads
\[
R_{k,m+1} = R_{k,m} + \frac{R_{k,m} - R_{k-1,m}}{4^{m} - 1}~.
\]
Implementation#
def romberg(
f,
a,
b,
accuracy=1e-8,
max_order=10
):
R = np.zeros((max_order, max_order))
h = (b - a) / 2.
R[0, 0] = h * (f(a) + f(b)) # The initial trapezoidal rule
for n in range(1, max_order):
trapezoid = 0.0
for j in range(2**(n-1)):
trapezoid += f(a + (2*j+1)*h)
R[n, 0] = 0.5 * R[n-1, 0] + h * trapezoid # The trapezoidal rule
l = 1
# The Romberg iterations
for m in range(1, n+1):
l *= 4
R[n, m] = (l * R[n, m-1] - R[n-1, m-1]) / (l-1)
print("Iteration: {0:5}, I = {1:20.15f}, error estimate = {2:10.15f}".format(n, R[n, m], abs(R[n, m] - R[n-1, m-1])))
if abs(R[n, m] - R[n-1, m-1]) < accuracy:
return R[n, m]
h /= 2.
print("Romberg method did not converge to required accuracy")
return R[-1, -1]
Illustration#
Let us compute our first integral using the Romberg method.
a = flimit_a
b = flimit_b
print("Computing the integral of",flabel, "over the interval (",a,",",b,") using Romberg method")
res = romberg(f,a,b,1e-6,18)
Computing the integral of x^4 - 2x + 2 over the interval ( 0.0 , 2.0 ) using Romberg method
Iteration: 1, I = 6.666666666666667, error estimate = 9.333333333333332
Iteration: 2, I = 6.400000000000000, error estimate = 0.266666666666667
Iteration: 3, I = 6.400000000000000, error estimate = 0.000000000000000
Now let us consider a more involved example of the Runge function
a = -2.
b = 2.
print("Computing the integral of",rungelabel, "over the interval (",a,",",b,") using Romberg method")
res =romberg(runge,a,b,1e-6,18)
Computing the integral of Runge function over the interval ( -2.0 , 2.0 ) using Romberg method
Iteration: 1, I = 2.679867986798680, error estimate = 2.640264026402640
Iteration: 2, I = 0.648895658796649, error estimate = 2.030972328002031
Iteration: 3, I = 0.554236075601252, error estimate = 0.094659583195396
Iteration: 4, I = 0.562270126297315, error estimate = 0.008034050696062
Iteration: 5, I = 0.587824850153293, error estimate = 0.025554723855978
Iteration: 6, I = 0.588636945021199, error estimate = 0.000812094867906
Iteration: 7, I = 0.588448788195693, error estimate = 0.000188156825505
Iteration: 8, I = 0.588451058525226, error estimate = 0.000002270329532
Iteration: 9, I = 0.588451069812733, error estimate = 0.000000011287507