Adaptive Simpson
CS/SE 4X03
Ned Nedialkov McMaster University February 24, 2021
Outline
Derivation of Simpson’s rule Adaptive Simpson Subtleties
Derivation of Simpson’s rule Adaptive Simpson Subtleties Derivation of Simpson’s rule
• Seek integration formula of the form
b a+b
f(x)dx ≈ Af(a) + Bf 2 + Cf(b) a
• Find A, B, C such that for quadratic polynomials the formula is exact:
b a+b
f(x)dx = Af(a) + Bf 2 + Cf(b) a
Copyright © 2021 N. Nedialkov. All rights reserved.
3/13
Derivation of Simpson’s rule Adaptive Simpson Subtleties Derivation of Simpson’s rule cont.
• Let a = −1, b = 1. We should integrate exactly 1, x, x2: 1
f(x) = 1 :
f(x) = x : f(x) = x2 :
dx = 2 = A + B + C xdx = 0 = −A + C
−1 12
x2dx = 3 = A + C
from which A = 1/3, C = 1/3, B = 4/3
• Hence
11
f(x)dx ≈ 3[f(−1) + 4f(0) + f(1)] −1
−1 1
−1
Copyright © 2021 N. Nedialkov. All rights reserved. 4/13
Derivation of Simpson’s rule Adaptive Simpson Subtleties Derivation of Simpson’s rule cont.
• Lety(x)=0.5(b−a)x+0.5(b+a),y(−1)=a,y(1)=b • Changing variables:
b b−a a+b
f(x)dx ≈ 6 f(a) + 4f 2 + f(b) a
Copyright © 2021 N. Nedialkov. All rights reserved.
5/13
Derivation of Simpson’s rule Adaptive Simpson Subtleties Adaptive Simpson
• Given a function f(x) on [a,b] and tolerance tol • find Q such that
where
|Q−I|≤ tol, b
I =
f(x)dx
a
Copyright © 2021 N. Nedialkov. All rights reserved. 6/13
Derivation of Simpson’s rule Adaptive Simpson Subtleties Adaptive Simpson cont.
Denoteh=b−a. Then b
where
I= f(x)dx=S(a,b)+E(a,b), a
h a + b S(a,b)=6 f(a)+4f 2 +f(b)
1 h5
E(a,b) = −90 2 f(4)(ξ), ξ between a and b
Denote S1 = S(a,b) and E1 = E(a,b)
Copyright © 2021 N. Nedialkov. All rights reserved. 7/13
Derivation of Simpson’s rule Adaptive Simpson Subtleties Adaptive Simpson cont.
• Letc=(a+b)/2andapplySimpsonon[a,c]and[c,b]: b
f(x)dx=S(a,c)+S(c,b)+E(a,c)+E(c,b)
S2 E2
• How to estimate the error? If f(4) does not change much on [a, b]
I =
• We can compute S1 and S2
a
1h/25 1 1h5
E(a,c) = −90 2 f(4)(ξ1) = 32 1 1h5
≈ − f(4)(ξ) 32 90 2
= 1 E1 32
−90 2
f(4)(ξ1)
,
ξ1 ∈ [a,c]
Copyright © 2021 N. Nedialkov. All rights reserved.
8/13
Derivation of Simpson’s rule Adaptive Simpson Subtleties Adaptive Simpson cont.
SimilarlyE(c,b)≈ 1E1 32
• Hence
E2=E(a,c)+E(c,b)≈ 1E1 16
• From I = S1 + E1 = S2 + E2,
S1 −S2 =E2 −E1 ≈E2 −16E2 =−15E2
• Then
E2≈ 1(S2−S1) 15
b1
f(x)dx = S2 + E2 ≈ S2 + 15(S2 − S1)
Copyright © 2021 N. Nedialkov. All rights reserved. 9/13
I =
a
Derivation of Simpson’s rule Adaptive Simpson Subtleties Method outline
Given f, [a, b] and tol :
• c=(a+b)/2
• ComputeS1 =S(a,b)andS2 =S(a,c)+S(c,b)
• E2 = (S2 − S1)/15
• If |E2| ≤ tol return S2 + E2
else apply recursively on [a, c] and [c, b]
Copyright © 2021 N. Nedialkov. All rights reserved.
10/13
Derivation of Simpson’s rule Adaptive Simpson Subtleties Adaptive Simpson cont.
Algorithm 2.1 (Adaptive Simpson). S = quadSimpson(f, a, b, tol)
h = b − a, c = (a + b)/2
S1 = h[f(a)+4f(a+b)+f(b)]
S2 = h [f(a)+4f(a+c)+2f(c)+4f(c+b)+f(b)] 12 2 2
E2 = 1 (S2 −S1) 15
if |E2| ≤ tol return S2 + E2
else
Q1 = quadSimpson(f, a, c, tol/2) Q2 = quadSimpson(f, c, b, tol/2) return Q1 + Q2
62
Copyright © 2021 N. Nedialkov. All rights reserved.
11/13
Derivation of Simpson’s rule Adaptive Simpson Subtleties Why it works
• LetI=bf(x)dx,I1 =cf(x)dx,I2 =cf(x)dx aab
• If|I1−Q1|≤ tol/2and|I2−Q2|≤ tol/2,then
|I − Q| = |(I1 + I2) − (Q1 + Q2)| = |I1 − Q1 + I2 − Q2|
≤ |I1 −Q1|+|I2 −Q2| ≤ tol/2+ tol/2
= tol
Copyright © 2021 N. Nedialkov. All rights reserved.
12/13
Derivation of Simpson’s rule Adaptive Simpson Subtleties Subtleties
• The error estimate assumes f(4) does not vary much, but it may, and then this estimate may not be accurate
To compensate in such cases, include a safety factor, e.g.
|E1| ≤ 10 × tol
• The recursion may run “deep” if tol is too small or f(4) varies a lot
Insert a counter to stop the recursion when the depth exceeds some number, e.g. 20
Copyright © 2021 N. Nedialkov. All rights reserved. 13/13