8.5 Peano kernel error formulas
#%config InlineBackend.figure_format = 'svg'
from pylab import *
The error formulae we have derived so far assume certain derivatives exist and are continuous. This may be a restrictive assumption in many cases. Even if the function is not differentiable everywhere, the derivatives may be integrable. In this case, we can use Taylor formula with an integral remainder term to perform the error analysis.
8.5.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 < ∞ Single interval. 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 ) \begin{align}
f(x) &= \clr{red}{f(a) + (x-a) f'(a)} + \clr{blue}{ \int_a^x (x-t) f''(t) \ud t } \\
&=: \clr{red}{p_1(x)} + \clr{blue}{ R_2(x) }
\end{align} 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 ( ∫ a b G ( x , t ) χ [ a , x ] ( t ) d t ) d x = ∫ a b ( ∫ a b G ( x , t ) χ [ a , x ] ( t ) d x ) d t = ∫ a b ( ∫ a b G ( x , t ) χ [ t , b ] ( x ) d x ) d t = ∫ a b ∫ t b G ( x , t ) d x d t \begin{align}
\int_a^b \int_a^x G(x,t) \ud t \ud x
&= \int_a^b \left( \int_a^b G(x,t) \chi_{[a,x]}(t) \ud t \right) \ud x \\
&= \int_a^b \left( \int_a^b G(x,t) \chi_{[a,x]}(t) \ud x \right) \ud t \\
&= \int_a^b \left( \int_a^b G(x,t) \chi_{[t,b]}(x) \ud x \right) \ud t \\
&= \int_a^b \int_t^b G(x,t) \ud x \ud t
\end{align} ∫ a b ∫ a x G ( x , t ) d t d x = ∫ a b ( ∫ a b G ( x , t ) χ [ a , x ] ( t ) d t ) d x = ∫ a b ( ∫ a b G ( x , t ) χ [ a , x ] ( t ) d x ) d t = ∫ a b ( ∫ a b G ( x , t ) χ [ t , b ] ( x ) d x ) d t = ∫ a b ∫ t b G ( x , t ) d x d t Here χ is the characteristic function and we used χ [ a , x ] ( t ) = χ [ t , b ] ( x ) \chi_{[a,x]}(t) = \chi_{[t,b]}(x) χ [ a , x ] ( t ) = χ [ t , b ] ( x ) a ≤ t ≤ x ≤ b a \le t \le x \le b a ≤ t ≤ x ≤ b G ( x , t ) = ( x − t ) f ′ ′ ( t ) G(x,t) = (x-t) f''(t) G ( x , t ) = ( x − t ) f ′′ ( t )
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 = ∫ a b ( t − a ) ( t − b ) 2 f ′ ′ ( t ) 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 \\
&= \int_a^b \frac{(t-a)(t-b)}{2} f''(t) \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 = ∫ a b 2 ( t − a ) ( t − b ) f ′′ ( t ) d t Composite rule. 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 if f ′ ′ f'' f ′′
1 b − a ∫ a b ∣ f ′ ′ ( t ) ∣ d t ≤ ∥ f ′ ′ ∥ ∞ \frac{1}{b-a} \int_a^b |f''(t)|\ud t \le \norm{f''}_\infty b − a 1 ∫ a b ∣ f ′′ ( t ) ∣ d t ≤ ∥ 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 ′ ′ ∥ ∞ = ∞ f''(x) = \frac{3}{4\sqrt{x}}, \qquad \norm{f''}_\infty = \infty f ′′ ( x ) = 4 x 3 , ∥ f ′′ ∥ ∞ = ∞ 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 )