Numerical Analysis Lecture (II): Numerical Integration

… views … visitors … total

中文版

You may want to read Numerical Analysis Lecture (I): Interpolation Methods first, because Newton-Cotes quadrature directly uses interpolation polynomials and Lagrange basis functions. This chapter then studies numerical integration: approximating definite integrals from finitely many function values.

---

Numerical Integration

This chapter discusses common methods for approximating definite integrals of the form

\[\int_a^b f(x)\,dx.\]

Integration task
Given an integrable function $f:[a,b]\to\mathbb{R}$, compute

\[I(f)=\int_a^b f(x)\,dx.\]

Even for functions that look relatively simple, an antiderivative may not be available in closed form, for example

\[\frac{\sin x}{x} \qquad\text{and}\qquad e^{-x^2}.\]

In such cases we need numerical integration. Many important quadrature methods are based on the same idea: choose several support points

\[(x_i,f(x_i)),\qquad x_i\in[a,b],\]

interpolate $f$ by a polynomial $p_n$, and then integrate that polynomial. This gives the approximation

\[I_n(f)=\int_a^b p_n(x)\,dx,\]

and methods of this type are called interpolatory quadrature.

2.1 Newton-Cotes Quadrature

2.1.1 Closed Newton-Cotes Quadrature

For $n\in\mathbb{N}$, choose equally spaced support points

\[x_i=a+ih,\qquad i=0,\ldots,n,\qquad h=\frac{b-a}{n}.\]

The Lagrange representation of the interpolation polynomial $p_n$ is

\[p_n(x)=\sum_{i=0}^{n} f(x_i)L_{i,n}(x), \qquad L_{i,n}(x)= \prod_{\substack{j=0\\ j\ne i}}^{n} \frac{x-x_j}{x_i-x_j}.\]

This gives the numerical quadrature formula

\[I_n(f) =\int_a^b p_n(x)\,dx =\sum_{i=0}^{n} f(x_i)\int_a^b L_{i,n}(x)\,dx.\]

With the substitution $x=a+sh$, where $s\in[0,n]$, we obtain:

Closed Newton-Cotes formula

\[I_n(f) =h\sum_{i=0}^{n}\alpha_{i,n}f(x_i), \qquad \alpha_{i,n} =\int_0^n \prod_{\substack{j=0\\ j\ne i}}^{n} \frac{s-j}{i-j}\,ds. \tag{2.1}\]

The numbers $\alpha_{0,n},\ldots,\alpha_{n,n}$ are called weights. They do not depend on $f$ or on the interval $[a,b]$, so they can be tabulated in advance. We always have

\[\sum_{i=0}^{n}h\alpha_{i,n}=b-a.\]

Definition 2.1.1
The integration formula

\[J(f)=\sum_{i=0}^{n}\beta_i f(x_i)\]

is said to have degree of precision $n$ if it integrates every polynomial of degree at most $n$ exactly.

By construction, the closed Newton-Cotes formula $I_n(f)$ has degree of precision $n$.

The next important question is how to estimate the error

\[E_n(f):=I(f)-I_n(f).\]

By Corollary 1.1.4,

\[|f(x)-p_n(x)| \le \frac{|f^{(n+1)}(\xi)|}{(n+1)!}(b-a)^{n+1},\]

where $\xi\in[a,b]$. Therefore, for some possibly different $\xi\in[a,b]$,

\[\left| \int_a^b f(x)\,dx-\int_a^b p_n(x)\,dx \right| \le \int_a^b |f(x)-p_n(x)|\,dx \le \frac{|f^{(n+1)}(\xi)|}{(n+1)!}(b-a)^{n+2}.\]

A Taylor expansion gives sharper remainder estimates. The results are listed below.

$n$	$\alpha_{i,n}$	Maximum error $E_n(f)$	Name
1	$\frac12,\frac12$	$-\frac{f^{(2)}(\xi)}{12}h^3$	Trapezoidal rule
2	$\frac13,\frac43,\frac13$	$-\frac{f^{(4)}(\xi)}{90}h^5$	Simpson’s rule
3	$\frac38,\frac98,\frac98,\frac38$	$-\frac{3f^{(4)}(\xi)}{80}h^5$	3/8 rule
4	$\frac{14}{45},\frac{64}{45},\frac{24}{45},\frac{64}{45},\frac{14}{45}$	$-\frac{8f^{(6)}(\xi)}{945}h^7$	Milne’s rule

One caveat is that for $n\ge 7$, Newton-Cotes weights can become negative. Negative weights introduce cancellation, which makes the formula increasingly unstable numerically. In practice, one usually does not try to improve accuracy simply by increasing $n$.

Diagram 2.1: Closed Newton-Cotes formulas include the endpoints. The trapezoidal rule connects the endpoint values by a line, while Simpson’s rule uses the quadratic curve through three nodes.

Closed Newton-Cotes: Trapezoidal Rule and Simpson's Rule

Closed formulas include $a$ and $b$ as nodes. A quadratic approximation often follows curvature better, but high-order Newton-Cotes formulas are not automatically more stable.

Example: Consider $f(x)=\log_2(x)$ on the interval $[a,b]=[2,4]$. The exact integral is

\[I(f) =\int_2^4 f(x)\,dx =\int_2^4 \log_2(x)\,dx =\left.\frac{1}{\ln(2)}(x\ln(x)-x)\right|_2^4\] \[=\frac{4\ln(4)-4-2\ln(2)+2}{\ln(2)} \approx 3.11461.\]

The trapezoidal rule gives $h=\frac{2}{1}=2$, hence

\[I_1(f) =2\left(\frac12 f(2)+\frac12 f(4)\right) =2\left(\frac12\log_2(2)+\frac12\log_2(4)\right) =2\left(\frac12+1\right)=3.\]

For Simpson’s rule, $h=\frac{2}{2}=1$, and

\[I_2(f) =1\left(\frac13 f(2)+\frac43 f(3)+\frac13 f(4)\right) \approx 3.11328.\]

For $\xi\in[2,4]$, the trapezoidal rule error estimate is

\[E_1(f) =-\frac{f''(\xi)h^3}{12} =\frac{1}{\xi^2\ln(2)}\frac{8}{12} \le \frac{1}{2^2\ln(2)}\frac{2}{3} \approx 0.24045.\]

This agrees reasonably well with the actual error $0.11461$. The Simpson error is

\[E_2(f) =-\frac{f^{(4)}(\xi)h^5}{90} =\frac{6}{\xi^4\ln(2)}\frac{1}{90} \le \frac{6}{2^4\ln(2)}\frac{1}{90} \approx 0.00601,\]

which agrees very well with the actual error $0.00133$.

2.1.2 Open Newton-Cotes Quadrature

For $n\in\mathbb{N}\cup{0}$, choose equally spaced support points inside the open interval $(a,b)$:

\[x_i=a+ih,\qquad i=1,\ldots,n+1,\qquad h=\frac{b-a}{n+2}.\]

Proceeding in the same way gives another interpolatory integration formula:

Open Newton-Cotes formula

\[\tilde I_n(f) =h\sum_{i=1}^{n+1}\tilde\alpha_{i,n}f(x_i), \qquad \tilde\alpha_{i,n} = \int_0^{n+2} \prod_{\substack{j=1\\ j\ne i}}^{n+1} \frac{s-j}{i-j}\,ds.\]

The quadrature error has a formula similar to the closed case:

$n$	$\tilde\alpha_{i,n}$	Maximum error $\tilde E_n(f)$	Name
0	$2$	$\frac{f^{(2)}(\xi)}{3}h^3$	Midpoint rectangle rule
1	$\frac32,\frac32$	$\frac{3f^{(2)}(\xi)}{4}h^3$
2	$\frac83,-\frac43,\frac83$	$\frac{28f^{(4)}(\xi)}{90}h^5$

Diagram 2.2: Open Newton-Cotes formulas do not use the endpoints. The midpoint rectangle rule samples only one interior point and multiplies its function value by the interval length.

Open Newton-Cotes: Midpoint Rectangle Rule

Open formulas are natural when endpoint values are unavailable, unreliable, or intentionally avoided by the algorithm.

Example: Again consider $f(x)=\log_2(x)$ on $[2,4]$, with exact value $I(f)\approx 3.11461$. The rectangle rule gives $h=\frac22=1$, hence

\[\tilde I_0(f)=2f(3)\approx 3.16992.\]

The error estimate is

\[\tilde E_0(f) =\frac{f^{(2)}(\xi)h^3}{3} =-\frac{1}{\xi^2\ln(2)}\frac13 \ge -\frac{1}{2^2\ln(2)}\frac13 \approx -0.12022,\]

while the exact error is $-0.05529$.

2.2 Composite Newton-Cotes Formulas

Newton-Cotes formulas are reliable mainly when the integration interval is small and the number of nodes is not too large. A more common strategy is therefore to split $[a,b]$ into smaller subintervals, apply a Newton-Cotes formula on each subinterval, and add up the resulting subintegrals.

Divide $[a,b]$ into $m$ subintervals of length

\[H=\frac{b-a}{m}.\]

Apply an $n$-degree Newton-Cotes formula on each subinterval and sum the results. Let

\[N=m\cdot n,\qquad H=\frac{b-a}{m},\qquad h=\frac{H}{n}=\frac{b-a}{N},\] \[x_i=a+ih,\qquad i=0,\ldots,N,\] \[y_j=a+jH,\qquad j=0,\ldots,m.\]

Since

\[I(f)=\sum_{j=0}^{m-1}\int_{y_j}^{y_{j+1}}f(x)\,dx,\]

we obtain:

Composite closed Newton-Cotes formula

\[S_N^{(n)}(f) =h\sum_{j=0}^{m-1}\sum_{i=0}^{n} \alpha_{i,n}f(x_{jn+i}).\]

The weights $\alpha_{i,n}$ are still given by (2.1). The quadrature error

\[R_N^{(n)}(f)=I(f)-S_N^{(n)}(f)\]

is obtained by adding the errors over all subintervals.

Theorem 2.2.1
Let $f\in C^{n+2}([a,b])$. Then there exists an intermediate point $\xi\in(a,b)$ such that

\[R_N^{(n)}(f)= \begin{cases} C(n)f^{(n+2)}(\xi)(b-a)h^{n+2}, & n \text{ is even},\\ C(n)f^{(n+1)}(\xi)(b-a)h^{n+1}, & n \text{ is odd}. \end{cases}\]

Here $C(n)$ is a constant depending only on $n$.

The most common composite formulas and their quadrature errors are listed below.

Composite trapezoidal rule

Closed, $n=1$, $h=\frac{b-a}{m}$:

\[S_N^{(1)}(f) =\frac{h}{2} \sum_{j=0}^{m-1} \left(f(x_j)+f(x_{j+1})\right), \qquad x_j=a+jh.\]

Error:

\[R_N^{(1)}(f) =-\frac{f''(\xi)}{12}(b-a)h^2.\]

Composite Simpson’s rule

Closed, $n=2$, $h=\frac{b-a}{2m}$:

\[S_N^{(2)}(f) =\frac{h}{3} \sum_{j=0}^{m-1} \left(f(x_{2j})+4f(x_{2j+1})+f(x_{2j+2})\right), \qquad x_j=a+jh.\]

Error:

\[R_N^{(2)}(f) =-\frac{f^{(4)}(\xi)}{180}(b-a)h^4.\]

Composite rectangle rule

Open, $n=0$, $2m=N$, $h=\frac{b-a}{N}$:

\[\tilde S_N^{(0)}(f) =2h\sum_{j=1}^{m} f(x_{2j-1}), \qquad x_j=a+jh.\]

Error:

\[\tilde R_N^{(0)}(f) =\frac{f''(\xi)}{6}(b-a)h^2.\]

Diagram 2.3: Composite formulas reduce the error by reducing the step size $h$. The leading error term of the trapezoidal rule usually scales like $h^2$, while Simpson’s rule scales like $h^4$.

Composite Formulas: Split the Interval and Sum Subintegrals

The composite idea is not to raise the interpolation degree without limit on one large interval, but to split the interval so that low-order formulas remain reliable on each small piece.

Example: Again consider $f(x)=\log_2(x)$ on $[2,4]$, where the exact value is $I(f)\approx 3.11461$, and take $m=2$.

• For the composite trapezoidal rule, $n=1$, $N=m=2$, and $h=\frac22=1$. Thus

  \[S_N^{(1)}(f) =\frac12\left(f(2)+f(3)+f(3)+f(4)\right) \approx \frac12(1+1.5849+1.5849+2) =3.08496.\]

  The error estimate gives

  \[R_N^{(1)}(f) =-\frac{f''(\xi)}{12}(b-a)h^2 = \frac{1}{\xi^2\ln(2)}\frac{2}{12} \le \frac{1}{2^2\ln(2)}\frac{2}{12} \approx 0.06011.\]

  This clearly reduces the error of the simple trapezoidal rule; the true error is $0.02965$.

• For the composite Simpson rule, $n=2$, $N=4$, and $h=\frac24=\frac12$. We get

  \[S_N^{(2)}(f) = \frac{1}{2\cdot 3} \left( (f(2)+4f(2.5)+f(3)) +(f(3)+4f(3.5)+f(4)) \right) \approx 3.11450.\]

  The error estimate is $R_N^{(2)}(f)\le 0.00037$, and the true error is $0.00011$.

• For the composite rectangle rule, $h=\frac24=\frac12$, and

  \[\tilde S_N^{(0)}(f) = \frac22\bigl(f(2.5)+f(3.5)\bigr) \approx 3.12928.\]

  The error estimate is

  \[\tilde R_N^{(0)}(f) = \frac{f''(\xi)}{6}(b-a)h^2 = -\frac{1}{\xi^2\ln(2)}\frac16\frac24 \ge -\frac{1}{2^2\ln(2)}\frac16\frac12 \approx -0.03005,\]

  and the true error is $-0.01467$.

Source, Copyright, and Usage Notes

This article mainly refers to the numerical analysis lecture notes in TU Darmstadt’s open repository: mathe3-script-2011-SoSe.pdf The upstream repository includes an Unlicense notice. This article is published for personal study, translation, and knowledge organization. The English wording, explanatory additions, and remade figures in this article do not represent the original authors or any official position. The personal organization, English text, explanatory notes, and remade figures in this article may be used for non-commercial study, discussion, and citation with attribution and the original link. Since part of this article is based on translation and organization of TU Darmstadt’s public lecture notes, the original material and any materials it may contain should remain subject to the original authors, repository, and license notices. For commercial use, systematic redistribution, publication, or large-scale adaptation, please verify the licensing status of the original material as well. If there are any translation, formula, terminology, or interpretation errors, or if the rights holder believes the material has been used improperly, please contact me and I will correct or remove it promptly.