Trapezoidal Rule - II - Numerical Analysis

from pylab import *
from scipy.integrate import fixed_quad

8.12.1Bernoulli polynomials¶

The Bernoulli polynomial $B_n(x)$ , $n \ge 0$ is defined implicitly by the generating function

\frac{t(\ee^{xt}-1)}{\ee^t - 1} = \sum_{j=1}^\infty B_j(x) \frac{t^j}{j!}

(1)

The first few polynomials are^[1]

B_0(x)=1, \qquad B_1(x)=x, \qquad B_2(x) = x^2 - x

(2)

B_3(x) = x^3 - \frac{3}{2} x^2 + \half x, \qquad B_4(x) = x^2(1-x)^2

(3)

Note that

B_k(0)=0, \qquad k \ge 1

(4)

It can be shown that for all $j > 1$

B_{2j}(x) > 0, \qquad 0 < x < 1

(5)

8.12.2Bernoulli numbers¶

These are defined implicitly by

\frac{t}{\ee^t - 1} = \sum_{j=0}^\infty B_j \frac{t^j}{j!}

(6)

The first few numbers are

B_0 = 1, \qquad B_1 = -\half, \qquad B_2 = \frac{1}{6}

(7)

B_4 = -\frac{1}{30}, \qquad B_6 = \frac{1}{42}, \qquad B_8 = -\frac{1}{30}

(8)

B_j = 0, \qquad \textrm{$j$ is odd and $j \ge 3$}

(9)

8.12.3Connection¶

To obtain the relation between $B_j(x)$ and $B_j$ , integrate the generating function of Bernoulli polynomial on $[0,1]$ to obtain

1 - \frac{t}{\ee^t - 1} = \sum_{j=1}^\infty \frac{t^j}{j!} \int_0^1 B_j(x) \ud x

(10)

Hence

B_0 = 1, \qquad B_j = -\int_0^1 B_j(x) \ud x, \quad j \ge 1

(11)

Next, define the periodic extension of $B_j(x)$ by

\bar{B}_j(x) = \begin{cases} B_j(x) & 0 \le x \le 1 \\ \bar{B}_j(x-1) & x > 1 \end{cases}

(12)

8.12.4Error of trapezoid rule¶

Theorem 1 (Error estimate)

Let $m \ge 0$ , $n \ge 1$ and define

h = (b-a)/n, \qquad x_j = a + jh, \qquad j=0,1,\ldots,n

(13)

Further assume that $f \in \cts^{2m+2}[a,b]$ for some $m \ge 0$ . Then the error in trapezoidal method is

\begin{aligned} E_n(f) &= \int_a^b f(x) \ud x - h {\sum_{j=0}^n}{}' f(x_j) \\ &= -\sum_{i=1}^m \frac{B_{2i} h^{2i}}{(2i)!}[ f^{(2i-1)}(b) - f^{(2i-1)}(a)] \\ & \quad + \frac{h^{2m+2}}{(2m+2)!} \int_a^b \bar{B}_{2m+2}((x-a)/h) f^{(2m+2)}(x) \ud x \end{aligned}

(14)

where the prime on the summation indicates that the first and last terms must be halved.

In more detail, the error of trapezoid method can be written as

\begin{align} E_n(f) &= \frac{h^2}{12} [ f^{(1)}(a) - f^{(1)}(b) ] - \frac{h^4}{720} [ f^{(3)}(a) - f^{(3)}(b) ] \\ & \qquad + \frac{h^6}{30240} [ f^{(5)}(a) - f^{(5)}(b) ] - \frac{h^8}{1209600} [ f^{(7)}(a) - f^{(7)}(b) ] + \ldots \end{align}

(15)

Proof 2

The periodicity of $f$ means

f^{(2j-1)}(a) = f^{(2j-1)}(b), \qquad j \ge 1

(16)

Hence for any $m$ the error formula of trapezoidal method yields

\begin{aligned} E_n(f) &= \frac{h^{2m+2}}{(2m+2)!} \int_a^b \bar{B}_{2m+2}((x-a)/h) f^{(2m+2)}(x) \ud x \\ &= - \frac{h^{2m+2}(b-a)B_{2m+2}}{(2m+2)!} f^{(2m+2)}(\xi), \qquad \xi \in [a,b] \\ &= \order{h^{2m+2}} \end{aligned}

(17)

where we used the positivity of $B_{2j}(x)$ and integral mean value theorem. Since $m$ was arbitrary, we obtain the desired result.

Example 1

Let us compute

\int_0^{2\pi} \exp(\sin x) \ud x

(18)

using trapezoid and Gauss quadrature.

f = lambda x: exp(sin(x))
a,b = 0.0,2*pi
qe = 7.954926521012844 # Exact integral

n,N = 2,10
e1,e2,nodes = zeros(N),zeros(N),zeros(N)
for i in range(N):
    x = linspace(a,b,n)
    val = trapezoid(f(x),dx=(b-a)/(n-1))
    e1[i] = abs(val - qe)
    val = fixed_quad(f,a,b,n=n)
    nodes[i] = n
    e2[i] = abs(val[0]-qe)
    print('%5d %20.10e %20.10e' % (n,e1[i],e2[i]))
    n = n+2

semilogy(nodes,e1,'o-')
semilogy(nodes,e2,'*-')
legend(('Trapezoid','Gauss'))
xlabel('n'), ylabel('Error');

    2     1.6717412138e+00     1.5277448438e+00
    4     2.8260084821e-04     2.0275305634e-01
    6     3.4594549447e-09     1.5148783594e-02
    8     7.1054273576e-15     1.2229354413e-03
   10     1.7763568394e-15     2.9323632517e-05
   12     1.7763568394e-15     3.4262054074e-06
   14     0.0000000000e+00     1.1795276755e-08
   16     2.6645352591e-15     4.7241393020e-09
   18     1.7763568394e-15     1.0296030695e-10
   20     1.7763568394e-15     2.8324009804e-12

The error of trapezoid becomes machine zero for 8 points while Gauss quadrature converges much more slowly.

Footnotes¶

$B_0(x)$ is defined to be one, it is not present in the definition of the generating function.
↩

References¶

Atkinson, K. E. (2004). An Introduction to Numerical Analysis (2nd ed.). Wiley.
Fornberg, B. (2021). Improving the Accuracy of the Trapezoidal Rule. SIAM Review, 63(1), 167–180. 10.1137/18M1229353

Numerical Analysis

Special techniques