Numerical Integration University of Massachusetts Amherst
The Trapezoidal Rule ¶ The simplest approach: approximate f ( x ) f(x) f ( x )
Single Interval ¶ ∫ a b f ( x ) d x ≈ b − a 2 ( f ( a ) + f ( b ) ) \int_a^b f(x)\,dx \approx \frac{b-a}{2}\left(f(a) + f(b)\right) ∫ a b f ( x ) d x ≈ 2 b − a ( f ( a ) + f ( b ) ) Connect ( a , f ( a ) ) (a, f(a)) ( a , f ( a )) ( b , f ( b ) ) (b, f(b)) ( b , f ( b ))
Local Error Analysis ¶ What is the error of this approximation? Let h = b − a h = b - a h = b − a
For f ∈ C 2 [ a , b ] f \in C^2[a, b] f ∈ C 2 [ a , b ] local truncation error (error on a single interval) is:
∫ a b f ( x ) d x − h 2 ( f ( a ) + f ( b ) ) = − h 3 12 f ′ ′ ( ξ ) \int_a^b f(x)\,dx - \frac{h}{2}\left(f(a) + f(b)\right) = -\frac{h^3}{12}f''(\xi) ∫ a b f ( x ) d x − 2 h ( f ( a ) + f ( b ) ) = − 12 h 3 f ′′ ( ξ ) for some ξ ∈ ( a , b ) \xi \in (a, b) ξ ∈ ( a , b ) h = b − a h = b - a h = b − a
Proof 1
Expand f ( x ) f(x) f ( x ) c = ( a + b ) / 2 c = (a+b)/2 c = ( a + b ) /2
f ( x ) = f ( c ) + f ′ ( c ) ( x − c ) + f ′ ′ ( η ( x ) ) 2 ( x − c ) 2 f(x) = f(c) + f'(c)(x - c) + \frac{f''(\eta(x))}{2}(x - c)^2 f ( x ) = f ( c ) + f ′ ( c ) ( x − c ) + 2 f ′′ ( η ( x )) ( x − c ) 2 Integrating from a a a b b b
∫ a b f ( x ) d x = h f ( c ) + f ′ ′ ( η ) 2 ∫ a b ( x − c ) 2 d x = h f ( c ) + h 3 24 f ′ ′ ( η ) \int_a^b f(x)\,dx = hf(c) + \frac{f''(\eta)}{2}\int_a^b (x-c)^2\,dx = hf(c) + \frac{h^3}{24}f''(\eta) ∫ a b f ( x ) d x = h f ( c ) + 2 f ′′ ( η ) ∫ a b ( x − c ) 2 d x = h f ( c ) + 24 h 3 f ′′ ( η ) where we used that ∫ a b ( x − c ) d x = 0 \int_a^b (x-c)\,dx = 0 ∫ a b ( x − c ) d x = 0 ∫ a b ( x − c ) 2 d x = h 3 / 12 \int_a^b (x-c)^2\,dx = h^3/12 ∫ a b ( x − c ) 2 d x = h 3 /12
For the trapezoidal approximation, expand f ( a ) f(a) f ( a ) f ( b ) f(b) f ( b ) c c c
f ( a ) = f ( c ) − h 2 f ′ ( c ) + h 2 8 f ′ ′ ( ξ 1 ) f(a) = f(c) - \frac{h}{2}f'(c) + \frac{h^2}{8}f''(\xi_1) f ( a ) = f ( c ) − 2 h f ′ ( c ) + 8 h 2 f ′′ ( ξ 1 ) f ( b ) = f ( c ) + h 2 f ′ ( c ) + h 2 8 f ′ ′ ( ξ 2 ) f(b) = f(c) + \frac{h}{2}f'(c) + \frac{h^2}{8}f''(\xi_2) f ( b ) = f ( c ) + 2 h f ′ ( c ) + 8 h 2 f ′′ ( ξ 2 ) Adding:
h 2 ( f ( a ) + f ( b ) ) = h f ( c ) + h 3 16 ⋅ f ′ ′ ( ξ 1 ) + f ′ ′ ( ξ 2 ) 2 \frac{h}{2}(f(a) + f(b)) = hf(c) + \frac{h^3}{16}\cdot\frac{f''(\xi_1) + f''(\xi_2)}{2} 2 h ( f ( a ) + f ( b )) = h f ( c ) + 16 h 3 ⋅ 2 f ′′ ( ξ 1 ) + f ′′ ( ξ 2 ) Taking the difference and using the intermediate value theorem to combine the f ′ ′ f'' f ′′
∫ a b f ( x ) d x − h 2 ( f ( a ) + f ( b ) ) = − h 3 12 f ′ ′ ( ξ ) \int_a^b f(x)\,dx - \frac{h}{2}(f(a) + f(b)) = -\frac{h^3}{12}f''(\xi) ∫ a b f ( x ) d x − 2 h ( f ( a ) + f ( b )) = − 12 h 3 f ′′ ( ξ ) for some ξ ∈ ( a , b ) \xi \in (a, b) ξ ∈ ( a , b )
Key observation: The local error is O ( h 3 ) O(h^3) O ( h 3 )
Composite Trapezoidal Rule ¶ For better accuracy, divide [ a , b ] [a, b] [ a , b ] n n n h = ( b − a ) / n h = (b-a)/n h = ( b − a ) / n x k = a + k h x_k = a + kh x k = a + kh k = 0 , 1 , … , n k = 0, 1, \ldots, n k = 0 , 1 , … , n
From Local to Global Error ¶ The global error is the total error when approximating the integral over [ a , b ] [a, b] [ a , b ]
For f ∈ C 2 [ a , b ] f \in C^2[a, b] f ∈ C 2 [ a , b ]
∫ a b f ( x ) d x − T n ( f ) = − ( b − a ) h 2 12 f ′ ′ ( ξ ) = O ( h 2 ) \int_a^b f(x)\,dx - T_n(f) = -\frac{(b-a)h^2}{12}f''(\xi) = O(h^2) ∫ a b f ( x ) d x − T n ( f ) = − 12 ( b − a ) h 2 f ′′ ( ξ ) = O ( h 2 ) for some ξ ∈ ( a , b ) \xi \in (a, b) ξ ∈ ( a , b )
Proof 2
On each subinterval [ x k − 1 , x k ] [x_{k-1}, x_k] [ x k − 1 , x k ]
∫ x k − 1 x k f ( x ) d x − h 2 ( f ( x k − 1 ) + f ( x k ) ) = − h 3 12 f ′ ′ ( ξ k ) \int_{x_{k-1}}^{x_k} f(x)\,dx - \frac{h}{2}(f(x_{k-1}) + f(x_k)) = -\frac{h^3}{12}f''(\xi_k) ∫ x k − 1 x k f ( x ) d x − 2 h ( f ( x k − 1 ) + f ( x k )) = − 12 h 3 f ′′ ( ξ k ) for some ξ k ∈ ( x k − 1 , x k ) \xi_k \in (x_{k-1}, x_k) ξ k ∈ ( x k − 1 , x k )
Summing over all n n n
∫ a b f ( x ) d x − T n ( f ) = − h 3 12 ∑ k = 1 n f ′ ′ ( ξ k ) \int_a^b f(x)\,dx - T_n(f) = -\frac{h^3}{12}\sum_{k=1}^{n} f''(\xi_k) ∫ a b f ( x ) d x − T n ( f ) = − 12 h 3 k = 1 ∑ n f ′′ ( ξ k ) Since f ′ ′ ∈ C [ a , b ] f'' \in C[a,b] f ′′ ∈ C [ a , b ] 1 n ∑ k = 1 n f ′ ′ ( ξ k ) \frac{1}{n}\sum_{k=1}^{n} f''(\xi_k) n 1 ∑ k = 1 n f ′′ ( ξ k ) min f ′ ′ \min f'' min f ′′ max f ′ ′ \max f'' max f ′′ ξ ∈ ( a , b ) \xi \in (a, b) ξ ∈ ( a , b )
1 n ∑ k = 1 n f ′ ′ ( ξ k ) = f ′ ′ ( ξ ) \frac{1}{n}\sum_{k=1}^{n} f''(\xi_k) = f''(\xi) n 1 k = 1 ∑ n f ′′ ( ξ k ) = f ′′ ( ξ ) Therefore:
∫ a b f ( x ) d x − T n ( f ) = − h 3 12 ⋅ n ⋅ f ′ ′ ( ξ ) = − h 3 n 12 f ′ ′ ( ξ ) \int_a^b f(x)\,dx - T_n(f) = -\frac{h^3}{12} \cdot n \cdot f''(\xi) = -\frac{h^3 n}{12}f''(\xi) ∫ a b f ( x ) d x − T n ( f ) = − 12 h 3 ⋅ n ⋅ f ′′ ( ξ ) = − 12 h 3 n f ′′ ( ξ ) Since n = ( b − a ) / h n = (b-a)/h n = ( b − a ) / h
∫ a b f ( x ) d x − T n ( f ) = − ( b − a ) h 2 12 f ′ ′ ( ξ ) \int_a^b f(x)\,dx - T_n(f) = -\frac{(b-a)h^2}{12}f''(\xi) ∫ a b f ( x ) d x − T n ( f ) = − 12 ( b − a ) h 2 f ′′ ( ξ ) Understanding Local vs Global Error ¶ Error Type Definition Trapezoidal Rule Local Error on one subinterval of width h h h O ( h 3 ) O(h^3) O ( h 3 ) Global Total error over [ a , b ] [a, b] [ a , b ] O ( h 2 ) O(h^2) O ( h 2 )
Why does the order drop from 3 to 2?
The global error accumulates local errors from n ∼ 1 / h n \sim 1/h n ∼ 1/ h
Global error ∼ n × Local error ∼ 1 h × h 3 = h 2 \text{Global error} \sim n \times \text{Local error} \sim \frac{1}{h} \times h^3 = h^2 Global error ∼ n × Local error ∼ h 1 × h 3 = h 2 This is the typical pattern: global order = local order − 1 .
Remark 1 (Python Implementation)
def trapezoidal(f, a, b, n):
"""Composite trapezoidal rule with n subintervals."""
h = (b - a) / n
x = np.linspace(a, b, n + 1)
y = f(x)
return h * (0.5 * y[0] + np.sum(y[1:-1]) + 0.5 * y[-1])Higher-Order Methods ¶ By using higher-degree polynomial approximations, we can achieve better accuracy.
Simpson’s rule uses a quadratic polynomial through three points ( a , f ( a ) ) (a, f(a)) ( a , f ( a )) ( m , f ( m ) ) (m, f(m)) ( m , f ( m )) ( b , f ( b ) ) (b, f(b)) ( b , f ( b )) m = ( a + b ) / 2 m = (a+b)/2 m = ( a + b ) /2
∫ a b f ( x ) d x ≈ h 6 ( f ( a ) + 4 f ( m ) + f ( b ) ) \int_a^b f(x)\,dx \approx \frac{h}{6}\left(f(a) + 4f(m) + f(b)\right) ∫ a b f ( x ) d x ≈ 6 h ( f ( a ) + 4 f ( m ) + f ( b ) ) where h = b − a h = b - a h = b − a O ( h 5 ) O(h^5) O ( h 5 ) O ( h 4 ) O(h^4) O ( h 4 )
Even higher-order Newton-Cotes formulas exist (using more equally-spaced points), though they become unstable for high orders. The optimal approach—Gaussian quadrature —chooses both the nodes and weights optimally and achieves remarkable efficiency. We will explore this in the chapter on interpolation .
Why Integration is Easier Than Differentiation ¶ Remark 2 (Smoothing vs. Roughening)
Operation Error behavior Differentiation Errors amplify (dividing by small h h h Integration Errors average out (summing many terms)
This is why numerical integration is generally more stable than numerical differentiation. Integration “smooths,” differentiation “roughens.”
From the perspective of conditioning:
This is why we can often integrate noisy data reliably, but differentiating noisy data is notoriously difficult.
Summary ¶ Rule Local Error Global Error Trapezoidal O ( h 3 ) O(h^3) O ( h 3 ) O ( h 2 ) O(h^2) O ( h 2 ) Simpson’s O ( h 5 ) O(h^5) O ( h 5 ) O ( h 4 ) O(h^4) O ( h 4 )
Key principle: Global order = local order − 1, because we sum O ( 1 / h ) O(1/h) O ( 1/ h )