Skip to main content
\(\require{mathrsfs}\require{cancel}\newcommand{\dee}[1]{\mathrm{d}#1} \newcommand{\half}{ \frac{1}{2} } \newcommand{\ds}{\displaystyle} \newcommand{\ts}{\textstyle} \newcommand{\es}{ {\varnothing}} \newcommand{\st}{ {\mbox{ s.t. }} } \newcommand{\pow}[1]{ \mathcal{P}\left(#1\right) } \newcommand{\set}[1]{ \left\{#1\right\} } \newcommand{\lin}{{\text{LIN}}} \newcommand{\quot}{{\text{QR}}} \newcommand{\simp}{{\text{SMP}}} \newcommand{\diff}[2]{ \frac{\mathrm{d}#1}{\mathrm{d}#2}} \newcommand{\bdiff}[2]{ \frac{\mathrm{d}}{\mathrm{d}#2} \left( #1 \right)} \newcommand{\ddiff}[3]{ \frac{\mathrm{d}^#1#2}{\mathrm{d}{#3}^#1}} \renewcommand{\neg}{ {\sim} } \newcommand{\limp}{ {\;\Rightarrow\;} } \newcommand{\nimp}{ {\;\not\Rightarrow\;} } \newcommand{\liff}{ {\;\Leftrightarrow\;} } \newcommand{\niff}{ {\;\not\Leftrightarrow\;} } \newcommand{\De}{\Delta} \newcommand{\bbbn}{\mathbb{N}} \newcommand{\bbbr}{\mathbb{R}} \newcommand{\bbbp}{\mathbb{P}} \newcommand{\cI}{\mathcal{I}} \newcommand{\cR}{\mathcal{R}} \newcommand{\cV}{\mathcal{V}} \newcommand{\Si}{\Sigma} \newcommand{\arccsc}{\mathop{\mathrm{arccsc}}} \newcommand{\arcsec}{\mathop{\mathrm{arcsec}}} \newcommand{\arccot}{\mathop{\mathrm{arccot}}} \newcommand{\erf}{\mathop{\mathrm{erf}}} \newcommand{\smsum}{\mathop{{\ts \sum}}} \newcommand{\atp}[2]{ \genfrac{}{}{0in}{}{#1}{#2} } \newcommand{\ave}{\mathrm{ave}} \newcommand{\llt}{\left \lt } \newcommand{\rgt}{\right \gt } \newcommand{\YEaxis}[2]{\draw[help lines] (-#1,0)--(#1,0) node[right]{$x$};\draw[help lines] (0,-#2)--(0,#2) node[above]{$y$};} \newcommand{\YEaaxis}[4]{\draw[help lines] (-#1,0)--(#2,0) node[right]{$x$};\draw[help lines] (0,-#3)--(0,#4) node[above]{$y$};} \newcommand{\YEtaxis}[4]{\draw[help lines] (-#1,0)--(#2,0) node[right]{$t$};\draw[help lines] (0,-#3)--(0,#4) node[above]{$y$};} \newcommand{\YEtaaxis}[4]{\draw[help lines, <->] (-#1,0)--(#2,0) node[right]{$t$}; \draw[help lines, <->] (0,-#3)--(0,#4) node[above]{$y$};} \newcommand{\YExcoord}[2]{\draw (#1,.2)--(#1,-.2) node[below]{$#2$};} \newcommand{\YEycoord}[2]{\draw (.2,#1)--(-.2,#1) node[left]{$#2$};} \newcommand{\YEnxcoord}[2]{\draw (#1,-.2)--(#1,.2) node[above]{$#2$};} \newcommand{\YEnycoord}[2]{\draw (-.2,#1)--(.2,#1) node[right]{$#2$};} \newcommand{\YEstickfig}[3]{ \draw (#1,#2) arc(-90:270:2mm); \draw (#1,#2)--(#1,#2-.5) (#1-.25,#2-.75)--(#1,#2-.5)--(#1+.25,#2-.75) (#1-.2,#2-.2)--(#1+.2,#2-.2);} \newcommand{\IBP}[7]{ \begin{array}{|c | l | l |} \hline \color{red}{\text{Option 1:}} & u=#2 &\color{red}{\dee{u}=#3 ~ \dee{#1}} \\ & \dee{v}=#5~\dee{#1} &\color{red}{v=#7} \\ \hline \color{blue}{\text{Option 2:}} & u=#5 &\color{blue}{\dee{u}=#6 ~ \dee{#1}} \\ &\dee{v}=#2 \dee{#1} &\color{blue}{v=#4} \\ \hline \end{array} } \renewcommand{\textcolor}[2]{{\color{#1}{#2}}} \newcommand{\trigtri}[4]{ \begin{tikzpicture} \draw (-.5,0)--(2,0)--(2,1.5)--cycle; \draw (1.8,0) |- (2,.2); \draw[double] (0,0) arc(0:30:.5cm); \draw (0,.2) node[right]{$#1$}; \draw (1,-.5) node{$#2$}; \draw (2,.75) node[right]{$#3$}; \draw (.6,1.1) node[rotate=30]{$#4$}; \end{tikzpicture}} \newcommand{\lt}{<} \newcommand{\gt}{>} \newcommand{\amp}{&} \)

Section1.11Numerical Integration

By now the reader will have come to appreciate that integration is generally quite a bit more difficult than differentiation. There are a great many simple-looking integrals, such as \(\int e^{-x^2}\dee{x}\text{,}\) that are either very difficult or even impossible to express in terms of standard functions  1 We apologise for being a little sloppy here — but we just want to say that it can be very hard or even impossible to write some integrals as some finite sized expression involving polynomials, exponentials, logarithms and trigonometric functions. We don't want to get into a discussion of computability, though that is a very interesting topic.. Such integrals are not merely mathematical curiosities, but arise very naturally in many contexts. For example, the error function

\begin{gather*} \mathrm{erf}(x) = \frac{2}{\sqrt{\pi}}\int_0^x e^{-t^2}\dee{t} \end{gather*}

is extremely important in many areas of mathematics, and also in many practical applications of statistics.

In such applications we need to be able to evaluate this integral (and many others) at a given numerical value of \(x\text{.}\) In this section we turn to the problem of how to find (approximate) numerical values for integrals, without having to evaluate them algebraically. To develop these methods we return to Riemann sums and our geometric interpretation of the definite integral as the signed area.

We start by describing (and applying) three simple algorithms for generating, numerically, approximate values for the definite integral \(\int_a^b f(x)\,\dee{x}\text{.}\) In each algorithm, we begin in much the same way as we approached Riemann sums.

  • We first select an integer \(n \gt 0\text{,}\) called the “number of steps”.
  • We then divide the interval of integration, \(a\le x\le b\text{,}\) into \(n\) equal subintervals, each of length \(\De x=\frac{b-a}{n}\text{.}\) The first subinterval runs from \(x_0=a\) to \(x_1=a+\De x\text{.}\) The second runs from \(x_1\) to \(x_2=a+2\De x\text{,}\) and so on. The last runs from \(x_{n-1}=b-\De x\) to \(x_n=b\text{.}\)

<<SVG image is unavailable, or your browser cannot render it>>

This splits the original integral into \(n\) pieces:

\begin{equation*} \int_a^b f(x)\,\dee{x} =\int_{x_0}^{x_1} f(x)\,\dee{x} +\int_{x_1}^{x_2} f(x)\,\dee{x} +\cdots +\int_{x_{n-1}}^{x_n} f(x)\,\dee{x} \end{equation*}

Each subintegral \(\int_{x_{j-1}}^{x_j} f(x)\,\dee{x}\) is approximated by the area of a simple geometric figure. The three algorithms we consider approximate the area by rectangles, trapezoids and parabolas (respectively).

<<SVG image is unavailable, or your browser cannot render it>>

We will explain these rules in detail below, but we give a brief overview here:

  1. The midpoint rule approximates each subintegral by the area of a rectangle of height given by the value of the function at the midpoint of the subinterval \begin{align*} \int_{x_{j-1}}^{x_{j}} f(x) \dee{x} & \approx f\left( \frac{x_{j-1}+x_{j}}{2} \right) \De x \end{align*} This is illustrated in the leftmost figure above.
  2. The trapezoidal rule approximates each subintegral by the area of a trapezoid with vertices at \((x_{j-1},0), (x_{j-1},f(x_{j-1})), (x_{j},f(x_{j})), (x_{j},0)\text{:}\) \begin{align*} \int_{x_{j-1}}^{x_{j}} f(x) \dee{x} & \approx \frac{1}{2} \left(f(x_{j-1})+f(x_j \right) \De x \end{align*} The trapezoid is illustrated in the middle figure above. We shall derive the formula for the area shortly.
  3. Simpson's rule approximates two adjacent subintegrals by the area under a parabola that passes through the points \((x_{j-1},f(x_{j-1}))\text{,}\) \((x_{j},f(x_{j}))\) and \((x_{j+1},f(x_{j+1}))\text{:}\) \begin{align*} \int_{x_{j-1}}^{x_{j+1}} f(x) \dee{x} & \approx \frac{1}{3} \left(f(x_{j-1})+4f(x_j)+f(x_{j+1} \right) \De x \end{align*} The parabola is illustrated in the right hand figure above. We shall derive the formula for the area shortly.

In what follows we need to refer to the midpoint between \(x_{j-1}\) and \(x_j\) very frequently. To save on writing (and typing) we introduce the notation

\begin{align*} \bar x_j &= \frac{1}{2} \left(x_{j-1}+x_j \right). \end{align*}