Numerical Integration – quantdev.blog

Interpolatory Quadrature Rules

Introduction

We are interested in an approximate calculation of the definite integral

\[ \begin{align*} I[f] = \int_{a}^{b}f(x)dx \end{align*} \tag{1}\]

where $f(x)$ is a given function and $[a,b]$ is a finite interval. This problem is often called numerical quadrature, since it relates to the ancient problem of the quadrature of a circle i.e. constructing a square with equal area to that of a circle. The computation of the above quantity is equivalent to solving the IVP:

\[ \begin{align*} y'(x) = f(x), \quad y(a)=0, \quad x \in[a,b] \end{align*} \tag{2}\]

for $y(b)=I[f]$.

As is well known, even many relatively simple integrals cannot be expressed in finite terms of elementary functions, and thus must be evaluated by numerical methods. Even when a closed form analytical solution exists, it may be preferable to use a numerical quadrature formula.

Since $I[f]$ is a linear functional, numerical integration is a special case of the problem of approximating a linear functional. We shall consider formulas of the form:

\[ \begin{align*} I[f] \approx \sum_{i=1}^n w_if(x_i) \end{align*} \tag{3}\]

where $x_1 < x_2 < \ldots < x_n$ are distinct nodes and $w_1$, $w_2$, $\ldots$, $w_n$ the corresponding weights. Often (but not always) all nodes lie in $[a,b]$.

Definition 1 (Order of accuracy of a Quadrature Rule) A quadrature rule (Equation 3) has order of accuracy (or degree of exactness) equal to $d$, iff it is exact for all polynomials of degree $\leq d$, that is, for all $p\in\mathcal{P}_{d+1}$.

Some Classical Formulas

Interpolatory quadrature formulas, where the nodes are constrained to be equally spaced, are called Newton-Cotes formulas. These are especially suited for integrating a tabulated function, a task that was more common before the computer age. The midpoint, the trapezoidal and the Simpson’s rules, to be described here, are all special cases of (unweighted) Newton-Cotes formulas.

The trapezoidal rule is based on the linear interpolation of $f(x)$ at $x_1 = a$ and $x_2 = b$, that is, $f(x)$ is approximated by :

\[ \begin{align*} p(x) = f(a) + (x-a)[a,b]f = f(a) + (x - a)\frac{f(b) - f(a)}{b - a} \end{align*} \]

The integral of $p(x)$ equals the area of a trapezoid with base $(b-a)$ times the average height $\frac{1}{2}(f(a) + f(b))$. Hence,

\[ \int_{a}^{b} f(x)dx \approx \frac{(b-a)}{2}(f(a) + f(b)) \]

To increase the accuracy, we subdivide the interval $[a,b]$ and assume that $f_i = f(x_i)$ is known on a grid of equidistant points:

\[ \begin{align*} x_0 = a, \quad x_i = x_0 + ih, \quad x_n = b \end{align*} \tag{4}\]

where $h = (b - a)/n$ is the step length. The trapezoidal approximation for the $i$th subinterval is:

\[ \begin{align*} \int_{x_i}^{x_{i+1}} f(x)dx = T(h) + R_i, \quad T(h) = \frac{h}{2}(f_i + f_{i+1}) \end{align*} \tag{5}\]

Let $p_2(x)\in\mathcal{P}_2$ be the unique interpolating polynomial (Newton polynomial) passing through the points $(x_i,f_i)$ and $(x_{i+1},f_{i+1})$, that is, $p_2(x_i)=f(x_i)$ and $p_2(x_{i+1}) = f(x_{i+1})$. The exact remainder in Newton’s interpolation formula is given by:

\[ \begin{align*} f(x) - p_2(x) &= [x_i,x_{i+1},x]f\cdot \Phi_2(x)\\ &=[x_i,x_{i+1},x]f \cdot (x - x_{i+1})\Phi_1(x)\\ &=[x_i,x_{i+1},x]f \cdot (x - x_{i+1})(x - x_{i})\Phi_0\\ &=[x_i,x_{i+1},x]f \cdot (x - x_{i+1})(x - x_{i}) \end{align*} \]

So, we have:

\[ \begin{align*} R_i &= \int_{x_i}^{x_{i+1}} (f(x) - p_2(x))dx = \int_{x_i}^{x_{i+1}}(x - x_i)(x - x_{i+1})[x_i,x_{i+1},x]f dx \end{align*} \tag{6}\]

By the theorem on the remainder term for interpolation, we can write:

\[ [x_1,\ldots,x_n,x_{n+1}]f = \frac{f^{(n)}(\xi)}{n!} \]

Consequently,

\[ [x_i,x_{i+1},x]f = \frac{f''(\xi)}{2} \]

So, we have:

\[ \begin{align*} R_i &= \frac{f''(\xi)}{2}\int_{x_i}^{x_{i+1}} (x - x_i)(x - x_{i+1})dx \end{align*} \tag{7}\]

Setting $x = x_i + ht$, $dx = hdt$ such that the limits of integration are from $t=0$ to $t=1$, $we get:

\[ \begin{align*} R_i &= \frac{f''(\xi)}{2} \int_{0}^{1}(ht)(x_i + ht - x_{i+1})h dt\\ &= \frac{f''(\xi)}{2} \int_{0}^{1}(ht)(ht - h)h dt\\ &= \frac{f''(\xi)}{2} \int_{0}^{1}(ht)(ht - h)h dt\\ &= \frac{f''(\xi) h^3}{2} \int_{0}^{1}(t^2 - t) dt\\ &= \frac{f''(\xi) h^3}{2} \left[\frac{t^3}{3} - \frac{t^2}{2}\right]_{0}^{1} \\ &= \frac{f''(\xi) h^3}{2} \left[\frac{1}{3} - \frac{1}{2}\right] \\ &= -\frac{1}{12}h^3 f''(\xi) \end{align*} \tag{8}\]

Summing the contributions for each subinterval $[x_i,x_{i+1}]$, $i=0...n$, gives:

\[ \begin{align*} \int_{a}^{b}f(x)dx = T(h) + R_T, \quad T(h) = \frac{h}{2}(f_0 + f_n) + h\sum_{i=1}^{n-1}f_i \end{align*} \tag{9}\]

which is the composite trapezoidal rule. The global truncation error is:

--- title: "Numerical Integration" author: "Quasar" date: "2024-11-11" categories: [Numerical Methods] image: "image.jpg" toc: true toc-depth: 3 format: html: code-tools: true code-block-border-left: true code-annotations: below highlight-style: pygments --- ## Interpolatory Quadrature Rules ### Introduction We are interested in an approximate calculation of the definite integral $$ \begin{align*} I[f] = \int_{a}^{b}f(x)dx \end{align*} $$ {#eq-definite-integral} where $f(x)$ is a given function and $[a,b]$ is a finite interval. This problem is often called **numerical quadrature**, since it relates to the ancient problem of the quadrature of a circle i.e. constructing a square with equal area to that of a circle. The computation of the above quantity is equivalent to solving the IVP: $$ \begin{align*} y'(x) = f(x), \quad y(a)=0, \quad x \in[a,b] \end{align*} $$ {#eq-IVP} for $y(b)=I[f]$. As is well known, even many relatively simple integrals cannot be expressed in finite terms of elementary functions, and thus must be evaluated by numerical methods. Even when a closed form analytical solution exists, it may be preferable to use a numerical quadrature formula. Since $I[f]$ is a linear functional, numerical integration is a special case of the problem of approximating a linear functional. We shall consider formulas of the form: $$ \begin{align*} I[f] \approx \sum_{i=1}^n w_if(x_i) \end{align*} $$ {#eq-quadrature-rules} where $x_1 < x_2 < \ldots < x_n$ are distinct **nodes** and $w_1$, $w_2$, $\ldots$, $w_n$ the **corresponding weights**. Often (but not always) all nodes lie in $[a,b]$. ::: {#def-quadrature-rule-of-order-of-accuracy-d} ### Order of accuracy of a Quadrature Rule A quadrature rule (@eq-quadrature-rules) has **order of accuracy** (or degree of exactness) equal to $d$, iff it is exact for all polynomials of degree $\leq d$, that is, for all $p\in\mathcal{P}_{d+1}$. ::: ### Some Classical Formulas Interpolatory quadrature formulas, where the nodes are constrained to be equally spaced, are called **Newton-Cotes** formulas. These are especially suited for integrating a tabulated function, a task that was more common before the computer age. The midpoint, the trapezoidal and the Simpson's rules, to be described here, are all special cases of (unweighted) Newton-Cotes formulas. The **trapezoidal** rule is based on the linear interpolation of $f(x)$ at $x_1 = a$ and $x_2 = b$, that is, $f(x)$ is approximated by : $$ \begin{align*} p(x) = f(a) + (x-a)[a,b]f = f(a) + (x - a)\frac{f(b) - f(a)}{b - a} \end{align*} $$ The integral of $p(x)$ equals the area of a trapezoid with base $(b-a)$ times the average height $\frac{1}{2}(f(a) + f(b))$. Hence, $$ \int_{a}^{b} f(x)dx \approx \frac{(b-a)}{2}(f(a) + f(b)) $$ To increase the accuracy, we subdivide the interval $[a,b]$ and assume that $f_i = f(x_i)$ is known on a grid of equidistant points: $$ \begin{align*} x_0 = a, \quad x_i = x_0 + ih, \quad x_n = b \end{align*} $$ {#eq-grid-points} where $h = (b - a)/n$ is the **step length**. The trapezoidal approximation for the $i$th subinterval is: $$ \begin{align*} \int_{x_i}^{x_{i+1}} f(x)dx = T(h) + R_i, \quad T(h) = \frac{h}{2}(f_i + f_{i+1}) \end{align*} $$ {#eq-trapezoidal-approx-in-ith-subinterval} Let $p_2(x)\in\mathcal{P}_2$ be the unique interpolating polynomial (Newton polynomial) passing through the points $(x_i,f_i)$ and $(x_{i+1},f_{i+1})$, that is, $p_2(x_i)=f(x_i)$ and $p_2(x_{i+1}) = f(x_{i+1})$. The exact remainder in Newton's interpolation formula is given by: $$ \begin{align*} f(x) - p_2(x) &= [x_i,x_{i+1},x]f\cdot \Phi_2(x)\\ &=[x_i,x_{i+1},x]f \cdot (x - x_{i+1})\Phi_1(x)\\ &=[x_i,x_{i+1},x]f \cdot (x - x_{i+1})(x - x_{i})\Phi_0\\ &=[x_i,x_{i+1},x]f \cdot (x - x_{i+1})(x - x_{i}) \end{align*} $$ So, we have: $$ \begin{align*} R_i &= \int_{x_i}^{x_{i+1}} (f(x) - p_2(x))dx = \int_{x_i}^{x_{i+1}}(x - x_i)(x - x_{i+1})[x_i,x_{i+1},x]f dx \end{align*} $$ {#eq-remainder-term} By the theorem on the remainder term for interpolation, we can write: $$ [x_1,\ldots,x_n,x_{n+1}]f = \frac{f^{(n)}(\xi)}{n!} $$ Consequently, $$ [x_i,x_{i+1},x]f = \frac{f''(\xi)}{2} $$ So, we have: $$ \begin{align*} R_i &= \frac{f''(\xi)}{2}\int_{x_i}^{x_{i+1}} (x - x_i)(x - x_{i+1})dx \end{align*} $$ {#eq-remainder-term-2} Setting $x = x_i + ht$, $dx = hdt$ such that the limits of integration are from $t=0$ to $t=1$, $we get: $$ \begin{align*} R_i &= \frac{f''(\xi)}{2} \int_{0}^{1}(ht)(x_i + ht - x_{i+1})h dt\\ &= \frac{f''(\xi)}{2} \int_{0}^{1}(ht)(ht - h)h dt\\ &= \frac{f''(\xi)}{2} \int_{0}^{1}(ht)(ht - h)h dt\\ &= \frac{f''(\xi) h^3}{2} \int_{0}^{1}(t^2 - t) dt\\ &= \frac{f''(\xi) h^3}{2} \left[\frac{t^3}{3} - \frac{t^2}{2}\right]_{0}^{1} \\ &= \frac{f''(\xi) h^3}{2} \left[\frac{1}{3} - \frac{1}{2}\right] \\ &= -\frac{1}{12}h^3 f''(\xi) \end{align*} $$ {#eq-error-in-subinterval} Summing the contributions for each subinterval $[x_i,x_{i+1}]$, $i=0...n$, gives: $$ \begin{align*} \int_{a}^{b}f(x)dx = T(h) + R_T, \quad T(h) = \frac{h}{2}(f_0 + f_n) + h\sum_{i=1}^{n-1}f_i \end{align*} $$ {#eq-composite-trapezoidal-rule} which is the **composite trapezoidal rule**. The global truncation error is: