Computational Fluid Dynamics - Elementary (i)

Integration: quadrature formulae

Integration is an important concept in mathematics and, together with its inverse, differentiation, is one of the two main operations in calculus. In this chapter, various integration techniques are explained and illustrated using Matlab codes.

Midpoint method

Applying this method it is possible to compute an approximation to an integral within an interval $ [a,b] $ by adding up the areas of a finite number of rectangles with heights defined by the function values, and user-defined widths which create sub-intervals (see Fig. i.1). To find the height of each rectangle, the function value at the midpoint of the sub-interval is used, which is more likely to be close to the average value of the function on the sub-interval. And the widths can be defined using linear, parabolic, cubic, etc, distributions. Then the midpoint method approximation to the integral for $ N $ rectangles can be given as:
\begin{equation} \int\limits_a^b f(x)dx \approx I_{md} = \sum_{i=1}^{N}f(\bar{x}_i)\Delta x_i, \end{equation} where:
\begin{equation} \bar{x}_i = \frac{x_{i-1}+x_i}{2} \ , \ \Delta x_i = x_i - x_{i-1}. \end{equation} For a linear distributed grid spacing, $ \Delta x = (b-a)/N $, the formula is reduced to:
\begin{equation} I_{md} = \frac{b-a}{N}\sum_{i=1}^{N}f(\bar{x}_i). \end{equation}

Fig. i.1: Midpoint approximation

Trapezoidal method

Assuming that the integrand $ f(x) $ is approximated as an $ n $-th order polynomial:
\begin{equation} f(x) \approx f_n(x) = a_0 + a_1 x + ... + a_{n-1} x^{n-1} + a_n x^n, \end{equation} and considering a linear polynomial $ n = 1 $, Trapezoidal approximation for a sub-interval $ [x_0,x_1] $ can be written as:
\begin{equation} \int\limits_{x_0}^{x_1} f(x)dx \approx \Delta x_1 \left[\frac{f(x_0) + f(x_1)}{2}\right]. \end{equation} For an interval $ [a,b] $ containing $ N $ sub-intervals we have (see Fig. i.2):
\begin{equation} \int\limits_{a}^{b} f(x)dx \approx I_{tr} = \sum_{i=1}^{N}\Delta x_i \left[\frac{f(x_{i-1}) + f(x_i)}{2}\right]. \end{equation}

Fig. i.2: Trapeziodal approximation

Simpson's method

In this method, instead of the linear segments used in the trapezoidal rule, parabolic arcs are used to approximate the integral of the function. Assuming the grid-points $ x_0 $, $ x_1 $ and $ x_2 $ in a equally spaced grid, the integral approximation over $ [x_0,x_2] $ is:
\begin{equation} \int\limits_{x_0}^{x_2} f(x)dx \approx \frac{\Delta x}{3} \left[f(x_0) + 4f(x_1) + f(x_2)\right]. \end{equation} Comparing Fig. i.3 with the previous approximation sketches shows that, for the given function, Simpson method results in the best approximation.

Fig. i.2: Trapeziodal approximation

The relevant Matlab codes can be downloaded using the following link:
Matlab Codes Bank