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.
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*} $$$$ L\left(f, p\right) \le \int_{a}^{b} f\left(x \right) \, \mathrm{d}x \le U\left(f, p\right). $$$$ \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*} $$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}$.
$$ \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)}$.
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)$.
$$ \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*} $$$$ R_i^1 = \dfrac{1}{3} \left( 4 R_i^0 - R_{i-1}^0 \right). $$$$ R_i^m = \dfrac{1}{4^m - 1} \left( 4^m R_i^{m-1} - R_{i-1}^{m-1} \right). $$Gaussian Quadrature
$$ \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$.
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}$.
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}$
$$ x = \dfrac{b-a}{2} t + \dfrac{a+b}{2}, $$$$ \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. $$