Numerical integration using trapezoidal and Simpson method

The derivation is provided for numerical approximation of the definite integral of a function $L(t)$ using trapezoidal and Simpson method over $t \in [t_0, ~t_f]$.

Time domain

Let $t \in [t_0, ~t_f]$ be divided into N equal intervals as follows, \([t_1=t_0,t_2,...,t_N,t_{N+1}=t_f]\) where, \(\Delta t_i=(t_{N+1}-t_1)/N=h~~ \forall i=1,...,N\)

Trapezoidal method

Normal

Approximate the integrand as a linear polynomial over $t \in [t_0, ~t_f]$

$\int_{t_0}^{t_f}L(t)dt=\dfrac{(t_f-t_0)}{2}(L(t_0)+L(t_f))$

However, this is not very accurate.

Composite

Approximate the integrand as a piecewise linear polynomial

$\int_{t_0}^{t_f}L(t)dt=\int_{t_1}^{t_2}L(t)dt+\int_{t_2}^{t_3}L(t)dt+…+\int_{t_{N-1}}^{t_N}L(t)dt+\int_{t_N}^{t_{N+1}}L(t)dt $

On $t \in [t_i ~t_{i+1}]$, $L$ is assumed to be linear and on each interval, the integral is the area of the trapezoid.

For $t \in [t_i~ t_{i+1}]$ and for $\tau \in [0~ h]$, the mappings between the domains is given by,

$\tau=0+(h-0)/(t_{i+1}+t_i)t=0+1t=t$

$d\tau=dt$

As they are equal, they can be used interchangeably. The $\tau$ domain simplifies the algebra and the linear approximation on the $i^{th}$ interval is given as

$L_i(\tau)=L_i+\dfrac{(L_{i+1}-L_i)}{\Delta \tau_i}\tau$

$\int L_i(\tau) d\tau=L_i\tau+ \dfrac{(L_{i+1}-L_i)}{\Delta \tau_i} \dfrac{\tau^2}{2}+C$

$\int_{0}^{h} L_i(\tau) d\tau=L_ih+ (L_{i+1}-L_i) \dfrac{h^2}{2h}=\dfrac{h}{2}(L_{i+1}+L_i), ~~~\Delta t_i=h $

This has to be summed over each interval to compute the intergal.

Simpson method

Normal

$\int_{t_0}^{t_f}L(t)dt=\dfrac{h}{6}\left(L(t_0)+4L(c=\dfrac{t_0+t_f}{2})+L(t_f)\right)$

The approximation of the integrand may not be accurate and it can lead to higher numerical errors in the integration.

Composite

Approximate integrand as a piecewise-quadratic polynomial. As the quadratic polynomial has 3 unknowns, they must be evaluated at 3 points. The $3^{rd}$ point is the midpoint of the interval denoted by $i+\frac{1}{2}$.

$L_i(\tau)=a\tau^2+b\tau+c$

$L_i=a0^2+b0+c=c$

$L_{i+1}=ah^2+bh+c=ah^2+bh+L_i$

$L_{i+1/2}=\dfrac{ah^2}{4}+\dfrac{bh}{2}+L_i$

In matrix form, the $2^{nd}$ and $3^{rd}$ equations can be written as, \(\begin{equation} \begin{bmatrix} a\\ b \end{bmatrix}=\dfrac{4}{(h^3)} \begin{bmatrix} \dfrac{h}{2}&-h\\ \dfrac{-h^2}{4}&{h}^{2} \end{bmatrix} \begin{bmatrix} L_{i+1}-L_i\\ L_{i+\frac{1}{2}}-L_i \end{bmatrix}= \begin{bmatrix} \dfrac{2}{h^2}&\dfrac{-4}{h^2}\\ \dfrac{-1}{h}&\dfrac{4}{h} \end{bmatrix} \begin{bmatrix} L_{i+1}-L_i\\ L_{i+\frac{1}{2}}-L_i \end{bmatrix} \end{equation}\)

The integral evaluates to as

$\int L_t(\tau)d\tau=\dfrac{a\tau^3}{3}+\dfrac{b\tau^2}{2}+c\tau+C$

$\int_{0}^{h} L_t(\tau)d\tau=\dfrac{ah^3}{3}+\dfrac{bh^2}{2}+ch$

\[\begin{aligned} \int_{0}^{h} L_t(\tau)d\tau &= L_ih+ \begin{bmatrix} \dfrac{h^3}{3} &\dfrac{h^2}{2} \end{bmatrix} \begin{bmatrix} a \\ b \end{bmatrix} \\ &=L_ih+ \begin{bmatrix} \dfrac{h^3}{3} &\dfrac{h^2}{2} \end{bmatrix} \begin{bmatrix} \dfrac{2}{h^2}&\dfrac{-4}{h^2}\\ \dfrac{-1}{h}&\dfrac{4}{h} \end{bmatrix} \begin{bmatrix} L_{i+1}-L_i\\ L_{i+\frac{1}{2}}-L_i \end{bmatrix}\\ &= L_ih+ h \begin{bmatrix} (\dfrac{2}{3}-\dfrac{1}{2}) & (2-\dfrac{4}{3}) \end{bmatrix} \begin{bmatrix} L_{i+1}-L_i\\ L_{i+\frac{1}{2}}-L_i \end{bmatrix}\\ &= L_ih+ h \begin{bmatrix} \dfrac{1}{6} & \dfrac{2}{3} \end{bmatrix} \begin{bmatrix} L_{i+1}-L_i\\ L_{i+\frac{1}{2}}-L_i \end{bmatrix}\\ &= h\left(\dfrac{6}{6}L_i+ -\dfrac{5}{6}L_i+\dfrac{1}{6}L_{i+1}+\dfrac{2}{3}L_{i+1}\right)\\ &= h\left(\dfrac{1}{6}L_i+\dfrac{1}{6}L_{i+1}+\dfrac{4}{6}L_{i+1/2}\right)\\ &=\dfrac{h}{6}\left(L_i+4L_{i+1/2}+L_{i+1}\right) \end{aligned}\]

Limitations

The method can be extended to higher degree polynomial but it will suffer from runge phenomemon and ill-conditioning. This can be avoided using optimally chosen non-uniform grid points.