Integration
Numerical Differentiation
Consider the approximations to the first-order derivative:
- Forward Difference Quotient: $$ \begin{equation*} D_{j}^{+} u = \dfrac{u_{j+1} - u_j}{h} \end{equation*} $$
- Backwards Difference Quotient: $$ \begin{equation*} D_{j}^{-} u = \dfrac{u_{j} - u_{j-1}}{h} \end{equation*} $$
- Central Difference Quotient: $$ \begin{equation*} D_{j}^{0} u = \dfrac{u_{j+1} - u_{j-1}}{2h} \end{equation*} $$ The forward and backwards difference schemes are first order approximations to the derivative. The central difference scheme is a second order accurate approximation.
Using the Taylor series of $f(x+h)$ and $f(x-h)$, the second-order derivative can be shown to be
$$ f^{\prime\prime} \left( x \right) = \dfrac{ f(x - h) - 2 f (x) + f(x+h)}{h^2} + \mathcal{O} \left( h^2 \right) . $$This is a method for deriving higher order approximations for derivatives from lower order approximations. Consider a first-order approximation to the derivative, $\varphi\left(h\right)$, such as backwards or forwards differencing, then
$$ f^{\prime}\left( x \right) = \varphi\left(h\right) + a_2 h^2 + a_3 h^3 + \ldots $$Now evaluate the derivative at ${h=h/2}$, so that
$$ f^{\prime}\left( x \right) = \varphi \left( h \right) + a_2 \left( \dfrac{h}{2} \right)^2 + a_3 \left( \dfrac{h}{2} \right)^3 + \ldots $$Combining the two terms so that the low order term cancel, i.e. via ${f^{\prime}\left( x \right) - 4 f^{\prime}\left( x \right)}$, then a better approximation can be found as
$$ f^{\prime}\left( x \right) = \varphi \left( h \right) - 4 \varphi \left( \dfrac{h}{2} \right) + \mathcal{O}\left(h^3\right). $$The process can also be applied to second order accurate schemes, such as central differencing, to produce more accurate approximations. The method can also be applied to higher order derivatives.
Numerical Integration
Create a partition, $p$, of the domain of integration: define $n+1$ nodes ${a = x_0 < x_1 < \ldots < x_n =b}$, so that there are $n$ sub-intervals $\left[x_i, x_{i+1} \right]$. Then approximate the area under the curve by summing the areas in each subinterval defined as
$$ \int\limits_{a}^{b}f\left( x\right) \, \mathrm{d}x \approx \sum\limits_{i=0}^{n-1} \left( x_{i+1} -x_i \right) f\left( x^{\ast} \right) \quad \textsf{where} \quad x^{\ast} \in \left[x_i, x_{i+1} \right]. $$If $f$ is continuous, the value of $x_{i}^{\ast}$ may be chosen arbitrarily in the interval $\left[ x_i, x_{i+1} \right]$. Then the lower and upper Riemann sums are given by
$$ \begin{align*} L\left(f, p\right) = & \sum\limits_{i=0}^{n-1} \left( x_{i+1} - x_i \right) m_i \quad \textsf{where} \quad m_i = \min_{x\in\left[ x_i, x_{i+1} \right]} f\left(x \right), \\ U\left(f, p\right) = & \sum\limits_{i=0}^{n-1} \left( x_{i+1} - x_i \right) M_i \quad \textsf{where} \quad M_i = \max_{x\in\left[ x_i, x_{i+1} \right]} f\left(x \right) \end{align*} $$so that
$$ L\left(f, p\right) \le \int_{a}^{b} f\left(x \right) \, \mathrm{d}x \le U\left(f, p\right). $$Additionally, the left and right Riemann sums are given, respectively, by
$$ \begin{align*} & \sum\limits_{i=0}^{n-1} \left( x_{i+1} -x_i \right) f\left(x_{i-1} \right), \\ & \sum\limits_{i=0}^{n-1} \left( x_{i+1} -x_i \right) f\left(x_{i} \right). \end{align*} $$Rather than rectangles, use trapezoids to approximate the integral in a sub domain
$$ \int\limits_{a}^{b}f\left( x\right) \mathrm{d}x \approx \sum\limits_{i=0}^{n-1} \left( x_{i+1} - x_i \right) \dfrac{ f\left( x_{i} \right) + f\left( x_{i+1} \right)}{2}. $$If the nodes of the partition are equally spaced, so that ${h=x_{i+1} - x_i}$, then the formula is given by
$$ T\left(f,p\right) = \dfrac{h}{2} \left( f\left(x_0 \right) + f\left(x_n \right) \right) + h \sum\limits_{i=1}^{n-2}f\left( x_{i} \right). $$Let ${f \in C^2\left( \left[a,b \right] \right)}$ and $p$ be equidistant partition of $\left[a,b \right]$, with ${h = x_{i+1} - x_i}$.
The error can be shown to have the form:
$$ \left| \int \limits_{a}^b f\left(x\right) \, \mathrm{d}x - T\left(f,p\right) \right| = a_2 h^2 + a_4 h^4 + \ldots $$that is, the error terms are even powers of the discretization. More precisely, the error for the trapezium rule can be written as
$$ \left| \int_{a}^b f\left(x\right) \, \mathrm{d}x - T\left(f,p\right) \right| = \dfrac{1}{12} \left| \left(b - a\right) h^2 f^{\prime\prime}\left( \xi \right) \right| $$for some ${\xi \in \left(a, b\right)}$.
The integral is approximated as
$$ \int_a^b f\left( x\right) \, \mathrm{d}x \approx \dfrac{b-a}{6} \left( f\left(a\right) + 4 f\left(\dfrac{a+b}{2}\right) +f\left(b\right) \right). $$Simpson’s rule uses interpolation with quadratic polynomials.
It can be applied in a composite manner, i.e. on many subdomains. It has an asymptotic error of $\mathcal{O}\left(h^4\right)$.
Let ${h=\dfrac{b-a}{2}}$, then Simpson’s Rule is written as
$$ \int_a^b f\left( x\right) \, \mathrm{d}x \approx \dfrac{h}{3} \left( f\left(a\right) + 4 f\left( a + h \right) +f\left(b\right) \right) $$and is referred to as Simpson’s 1/3 Rule. It uses quadratic interpolation of the function through the end and midpoints.
Romberg’s method uses the Trapezoidal Rule and then Richardson Extrapolation to estimate integrals. First consider a sequence of partitions, $p_i$, for ${i=0, \ldots, n}$, of equal spacing given by ${h_i = \dfrac{b-a}{2^i}}$ which yield a sequence of integrals ${R_i^0 = T_i \left( f, p_i \right)}$. Refinements of the integrals can then be produced by Richardson Extrapolation:
$$ \begin{array}{cc} R_0^0 & \\ & R_1^1 \\ R_1^0 & \\ & R_2^1 \\ R_2^0 & \\ \vdots & \\ R_n^0 & \end{array} $$To evaluate the next set of values, consider the two integrals
$$ \begin{align*} \int_{a}^{b} f\left( x \right) \, \mathrm{d} x & = R_{i-1}^0 + a_2 h^2 + a_4 h^4 + \ldots \\ \int_{a}^{b} f\left( x \right) \, \mathrm{d} x & = R_{i}^0 + a_2 \left( \dfrac{h}{2} \right)^2 + a_4 \left( \dfrac{h}{2} \right)^4 + \ldots \end{align*} $$Note that there no odd terms in the error. Then, define the next approximation as
$$ R_i^1 = \dfrac{1}{3} \left( 4 R_i^0 - R_{i-1}^0 \right). $$which has an error $\mathcal{O}\left(h^4\right)$. The extrapolated values are equivalent to integrals approximated by Simpson’s rule. The recurrence formula can be derived
$$ R_i^m = \dfrac{1}{4^m - 1} \left( 4^m R_i^{m-1} - R_{i-1}^{m-1} \right). $$Gaussian Quadrature
Generalise the quadrature formula so that an integral is approximated as
$$ \begin{equation*} I_{n}[f] = \sum \limits_{i=0}^{n} A_{i} f\left( x_{i} \right). \end{equation*} $$The above equation is a weighted sum of the values of $f$ at the points $x_{i}$, for $i=0, \ldots, n$. These points are said to be the nodes of the quadrature formula, while the $A_{i} \in \mathbb{R}$ are its coefficients or weights. Both weights and nodes depend in general on $n$.
- Can the weights be chosen such that the error in an integral is minimized?
- Furthermore, can the nodes be chosen such that the integral can be further improved?
Two real-valued functions $f(x)$ and $g(x)$ are said to be orthogonal if
$$ \langle f, g \rangle = \int_a^b f(x) g(x) \mathrm{d} x = 0. $$where $\langle f, g \rangle$ satisfies the conditions to be an inner product.
Let $q\left(x\right)$ be a non-trivial polynomial of degree ${n+1}$ such that
- it has $n+1$ distinct roots, denoted as $x_i$, in $\left[ a, b \right]$,
- the polynomial $q(x)$ satisfies $$ \int_{a}^{b} x^k q\left(x\right) \, \mathrm{d}x = 0 \quad \textsf{for} \quad k=0,\ldots, n. $$ i.e. $q(x)$ is orthogonal to $x^k$.
Then, denote the integral as
$$ I_n\left[ f\right] = \int_{a}^{b} f\left(x\right) \, \mathrm{d}x = \sum\limits_{i=0}^n A_i f\left( x_i \right) $$with nodes given by the roots of the orthogonal polynomial $q$, and weights given by $A_i = \int_{a}^{b} l_i \left(x\right) \, \mathrm{d}x$, using Lagrange interpolation for all polynomials $f\left(x\right)$ of degree less than or equal to ${2n+1}$.
Then, the integral $I_n\left[ f\right]$ integrates all polynomials of degree ${2n+1}$ exactly.
It is said that the degree of exactness of $I\left[ f\right]$ is ${2n+1}$.
The Legendre polynomials are a set of orthogonal polynomials where
$$ \int_{-1}^1 P_m(x) P_n(x) \, \mathrm{d} x = 0 \quad \textsf{for} \quad n \ne m. $$and $P_0 = 1$. Thus, $P_1 = x$, ${P_2 = \left(3x^2 -1\right) / 2}$. The Legendre polynomials obey a recursive formula:
$$ P_n = \dfrac{2n-1}{n} x P_{n-1}(x) - \dfrac{n-1}{n} P_{n-2}(x) \quad \textsf{for} \quad n \ge 2. $$Gauss-Legendre quadrature uses the roots of the Legendre polynomials as the nodes for integration, and weights may also be found by equating the quadrature expressions with the exact integrals for orders of $f(x)$, i.e. ${f(x)=1, x, x^2, \ldots}$
The domain of integration can be scaled from $(-1,1)$ to $(a, b)$ via the invertible transformation
$$ x = \dfrac{b-a}{2} t + \dfrac{a+b}{2}, $$so that
$$ \int_a^b f(x) \, \mathrm{d}x = \dfrac{b - a}{2} \int_{-1}^1 f\left( \dfrac{b-a}{2} t + \dfrac{a+b}{2}\right) \, \mathrm{d}t. $$