The error formulae we have derived so far assume certain derivatives exist and are continuous. This may be a restrictive assumption in some cases. Even if the function is not differentiable everywhere, it the derivatives may be integrable. In this case, we can Taylor formula with an integral remainder term to perform the error analysis.
1 Trapezoidal rule ¶ Assume that f f f f ′ ′ f'' f ′′ [ a , b ] [a,b] [ a , b ]
∫ a b ∣ f ′ ′ ( x ) ∣ d x < ∞ \int_a^b |f''(x)| \ud x < \infty ∫ a b ∣ f ′′ ( x ) ∣ d x < ∞ Then Taylor’s theorem gives
f ( x ) = f ( a ) + ( x − a ) f ′ ( a ) + ∫ a x ( x − t ) f ′ ′ ( t ) d t = : p 1 ( x ) + R 2 ( x ) f(x) = f(a) + (x-a) f'(a) + \int_a^x (x-t) f''(t) \ud t =: p_1(x) + R_2(x) f ( x ) = f ( a ) + ( x − a ) f ′ ( a ) + ∫ a x ( x − t ) f ′′ ( t ) d t =: p 1 ( x ) + R 2 ( x ) The quadrature error E n E_n E n
E n ( F + G ) = E n ( F ) + E n ( G ) , F , G ∈ C [ a , b ] E_n(F+G) = E_n(F) + E_n(G), \qquad F,G \in \cts[a,b] E n ( F + G ) = E n ( F ) + E n ( G ) , F , G ∈ C [ a , b ] Thus
E 1 ( f ) = E 1 ( p 1 ) + E 1 ( R 2 ) = E 1 ( R 2 ) E_1(f) = E_1(p_1) + E_1(R_2) = E_1(R_2) E 1 ( f ) = E 1 ( p 1 ) + E 1 ( R 2 ) = E 1 ( R 2 ) since Trapezoidal rule is exact for linear polynomials. Now
E 1 ( R 2 ) = ∫ a b R 2 ( x ) d x − b − a 2 [ R 2 ( a ) + R 2 ( b ) ] , R 2 ( a ) = 0 = ∫ a b ∫ a x ( x − t ) f ′ ′ ( t ) d t d x − b − a 2 ∫ a b ( b − t ) f ′ ′ ( t ) d t \begin{aligned}
E_1(R_2)
&= \int_a^b R_2(x) \ud x - \frac{b-a}{2}[R_2(a) + R_2(b)], \qquad R_2(a) = 0 \\
&= \int_a^b \int_a^x (x-t) f''(t) \ud t \ud x - \frac{b-a}{2}\int_a^b (b-t) f''(t) \ud t
\end{aligned} E 1 ( R 2 ) = ∫ a b R 2 ( x ) d x − 2 b − a [ R 2 ( a ) + R 2 ( b )] , R 2 ( a ) = 0 = ∫ a b ∫ a x ( x − t ) f ′′ ( t ) d t d x − 2 b − a ∫ a b ( b − t ) f ′′ ( t ) d t For any integrable function G ( x , t ) G(x,t) G ( x , t )
∫ a b ∫ a x G ( x , t ) d t d x = ∫ a b ∫ t b G ( x , t ) d x d t \int_a^b \int_a^x G(x,t) \ud t \ud x = \int_a^b \int_t^b G(x,t) \ud x \ud t ∫ a b ∫ a x G ( x , t ) d t d x = ∫ a b ∫ t b G ( x , t ) d x d t Using this result in the above error formula
E 1 ( R 2 ) = ∫ a b f ′ ′ ( t ) ∫ t b ( x − t ) d x d t − b − a 2 ∫ a b ( b − t ) f ′ ′ ( t ) d t = 1 2 ∫ a b f ′ ′ ( t ) ( t − a ) ( t − b ) d t \begin{aligned}
E_1(R_2)
&= \int_a^b f''(t) \int_t^b (x-t) \ud x \ud t - \frac{b-a}{2} \int_a^b (b-t) f''(t) \ud t \\
&= \half \int_a^b f''(t) (t-a)(t-b) \ud t
\end{aligned} E 1 ( R 2 ) = ∫ a b f ′′ ( t ) ∫ t b ( x − t ) d x d t − 2 b − a ∫ a b ( b − t ) f ′′ ( t ) d t = 2 1 ∫ a b f ′′ ( t ) ( t − a ) ( t − b ) d t For the composite Trapezoid rule, the error is obtained by adding the error from each interval and can be written as
E n ( f ) = ∫ a b K ( t ) f ′ ′ ( t ) d t E_n(f) = \int_a^b K(t) f''(t) \ud t E n ( f ) = ∫ a b K ( t ) f ′′ ( t ) d t where
K ( t ) = 1 2 ( t − x j − 1 ) ( t − x j ) , t ∈ [ x j − 1 , x j ] , j = 1 , 2 , … , n K(t) = \half (t - x_{j-1})(t - x_j), \qquad t \in [x_{j-1}, x_j], \qquad j=1,2,\ldots,n K ( t ) = 2 1 ( t − x j − 1 ) ( t − x j ) , t ∈ [ x j − 1 , x j ] , j = 1 , 2 , … , n This gives the following error estimate
∣ E n ( f ) ∣ ≤ ∥ K ∥ ∞ ∫ a b ∣ f ′ ′ ( t ) ∣ d t = h 2 8 ∫ a b ∣ f ′ ′ ( t ) ∣ d t |E_n(f)| \le \norm{K}_\infty \int_a^b |f''(t)| \ud t = \frac{h^2}{8} \int_a^b |f''(t)| \ud t ∣ E n ( f ) ∣ ≤ ∥ K ∥ ∞ ∫ a b ∣ f ′′ ( t ) ∣ d t = 8 h 2 ∫ a b ∣ f ′′ ( t ) ∣ d t Compare this to the previous error estimate
∣ E n ( f ) ∣ = b − a 12 h 2 ∣ f ′ ′ ( η ) ∣ ≤ b − a 12 h 2 ∥ f ′ ′ ∥ ∞ |E_n(f)| = \frac{b-a}{12} h^2 |f''(\eta)| \le \frac{b-a}{12} h^2 \norm{f''}_\infty ∣ E n ( f ) ∣ = 12 b − a h 2 ∣ f ′′ ( η ) ∣ ≤ 12 b − a h 2 ∥ f ′′ ∥ ∞ which may over-estimate the error[1] f ′ ′ f'' f ′′
The function
f ( x ) = x 3 / 2 , x ∈ [ 0 , 1 ] f(x) = x^{3/2}, \qquad x \in [0,1] f ( x ) = x 3/2 , x ∈ [ 0 , 1 ] does not have finite second derivative since
f ′ ′ ( x ) = 3 4 x f''(x) = \frac{3}{4\sqrt{x}} f ′′ ( x ) = 4 x 3 The error formula based on derivative cannot be used but since
∫ 0 1 ∣ f ′ ′ ( x ) ∣ d x = 3 2 \int_0^1 |f''(x)|\ud x = \frac{3}{2} ∫ 0 1 ∣ f ′′ ( x ) ∣ d x = 2 3 the error formula based on integral remainder term can be used to conclude that
E n ( f ) ≤ 3 h 2 16 → 0 as n → ∞ E_n(f) \le \frac{3 h^2}{16} \to 0 \qquad \textrm{as $n \to \infty$} E n ( f ) ≤ 16 3 h 2 → 0 as n → ∞ # Performs Trapezoid quadrature
def integrate(a,b,n,f):
h = (b-a)/n
x = linspace(a,b,n+1)
y = f(x)
res1 = 0.5*y[0] + sum(y[1:n]) + 0.5*y[n]
return h*res1
f = lambda x: x * sqrt(x)
a, b = 0.0, 1.0
Ie = 2.0/5.0 # Exact integral
n, N = 2, 10
e = zeros(N)
for i in range(N):
h = (b-a)/n
I = integrate(a,b,n,f)
e[i] = Ie - I
if i > 0:
print('%6d %18.8e %14.5g %18.8e'% (n,e[i],e[i-1]/e[i],3*h**2/16))
else:
print('%6d %18.8e %14.5g %18.8e'%(n,e[i],0,3*h**2/16))
n = 2*n 2 -2.67766953e-02 0 4.68750000e-02
4 -7.01811086e-03 3.8154 1.17187500e-02
8 -1.81246480e-03 3.8721 2.92968750e-03
16 -4.63401302e-04 3.9112 7.32421875e-04
32 -1.17671210e-04 3.9381 1.83105469e-04
64 -2.97398625e-05 3.9567 4.57763672e-05
128 -7.49190899e-06 3.9696 1.14440918e-05
256 -1.88304417e-06 3.9786 2.86102295e-06
512 -4.72540682e-07 3.9849 7.15255737e-07
1024 -1.18449772e-07 3.9894 1.78813934e-07
We still get O ( h 2 ) O(h^2) O ( h 2 )