Numerical integration

When to use numerical integration?

  • Non-integrable Functions: In many real-world applications, you'll encounter functions that are not easily integrable in closed-form, making numerical methods the go-to strategy.
  • Computational Efficiency: Even when an analytical solution exists, its complexity might render it computationally inefficient. In such cases, a numerical approximation can often provide sufficient accuracy more quickly.
  • Data-Driven Models: When dealing with empirical data that doesn't conform to an analytical model, numerical integration is an effective way to derive value from such data.
    • Electrical Systems: For instance, you may need to integrate the power consumption curve of a machine over time to calculate total energy usage.
    • Computational Geometry: When you have a complex object represented as a 3D mesh, numerical integration can help in calculating its volume, surface area, or other geometric properties.

Method

Start from the mathematical definition and split into $N$ intervals

$$ \int_a^b f(x) dx = \sum_{i=0}^{N-1} \int\limits_{x_i}^{x_{i+1}} f(x) dx $$

we have $x_0=a$, $x_N=b$.

We define $\Delta x = \frac{b-a}{N}$.

If panels are narrow enough we can approximate the function inside the panel $i$:

$$ f(x) = f(x_i) + f'(x_i) h + \frac12f''(x_i) h^2 +... $$

with $h = x-x_i$.

We can now devise different integration methods with differring degree of precision by including more of fewer terms in the expansion.

Rectangle method

For the rectangle method we only take the first term:

$$ \int_a^b f(x) dx = \sum_{i=0}^{N-1} \int\limits_{x_i}^{x_{i+1}} f(x_i)dx \approx \frac{b-a}{N}\sum_{i=0}^{N-1}f(x_i) $$

We can also expand the function from the right boundary of each panel and we get

$$ \int_a^b f(x) dx \approx \sum_{i=0}^{N-1} \int\limits_{x_i}^{x_{i+1}} f(x_i) dx \approx \frac{b-a}{N}\sum_{i=1}^{N}f(x_i) $$$$= \Delta x \sum\limits_{i=1}^{N}f(x_i)$$
right rectangles left rectangles

Or we can expand around the middle of the interval:

$$ \int_a^b f(x) dx = \sum_{i=0}^{N-1} \int\limits_{x_i}^{x_{i+1}} f(m_i) dx \approx \frac{b-a}{N}\sum_{i=1}^{N}f(m_i) = \Delta x \sum_{i=1}^{N}f(m_i) $$

where

$$m_i = \frac{x_i+x_{i-1}}{2}$$

Error estimate

In "normal" conditions the error scales with increasing number of panels $N$ as

  • $N^{-1}$ for the left- or right-rectangle method
  • $N^{-2}$ for the midpoint method

Trapezium rule

We can improve the situation by increasing the order of the approximation. If we include the first derivative information we obtain the trapezium rule, which approximates the curve with linear segments:

Derivation

$$ \int_a^b f(x) dx \approx \sum_{i=0}^{N-1} \int\limits_{x_i}^{x_{i+1}} f(x_i)+ (x-x_i)\frac{f(x_{i+1})-f(x_i)}{x_{i+1}-x_1} $$$$ = \sum_{i=0}^{N-1} \left( \Delta x f(x_i) + \frac{1}{2}\Delta x^2\frac{f(x_{i+1})-f(x_i)}{\Delta x}\right) $$$$ = \sum_{i=0}^{N-1} \Delta x \frac12 \left( f(x_i)+f(x_{i+1}) \right) =\Delta x\left[\sum_{i=1}^{N-1}f(x_i) +\frac12\left(f(a)+f(b)\right)\right]$$

Cost

  • requires $N+1$ function evaluation
  • one more than the rectangle method
  • error scales like $N^{-2}$

Error estimate

Simpson Rule

  • For each panel we can approximate the function by a quadratic function.
  • we need three points for each panel

Derivation

Ansatz:

$$g(x)=a x^2 + bx +c$$

Solve equations

$ g(0) = f_1 $,

$ g(h) = f_2 $,

$ g(2h) = f_3 $ for $a$, $b$, $c$

Integrate $g(x)$ over the interval $[0,2h]$ gives:

$$ \frac{h}{3}\left(f_1+4f_2+f_3\right) = \Delta x\frac{1}{6}\left(f_1+4f_2+f_3\right). $$
In [1]:
import sympy
x, h, f1, f2, f3 = sympy.symbols("x h f1 f2 f3")
a, b, c = sympy.symbols("a b c")

def g(x):
    return a*x**2 + b*x + c

abcRule = sympy.solve([
    g(0) - f1,
    g(h) - f2,
    g(2 * h) - f3,
],a,b,c)

integ = sympy.integrate(g(x) , (x , 0, 2*h ) )

integ.subs(abcRule).simplify()
Out[1]:
$\displaystyle \frac{h \left(f_{1} + 4 f_{2} + f_{3}\right)}{3}$

Now combine all panels:

$$ \int_a^b f(x) dx \approx \frac{ \Delta x}{6} \bigg( f(x_0) + 4 f(m_1) + 2 f(x_1) + 4 f(m_2) + 2 f(x_2) + ... + 4 f(m_{N-1})+ 2 f(x_{N-1})+ 4 f(m_{N})+ f(x_N) \bigg)$$

Error estimate

Higher order approximation

  • we can use higher and higher order approximations
  • faster convergence for smooth functions
  • numerical round-off issues

Validity

  • Higher order only helps if the function and its derivatives are smooth
$$f(x) = |x - \pi |\qquad I=\int\limits_0^6 f(x)\,dx$$

  • Only help if there are no divergences
$$ f(x) = \sqrt{1-x^2},\qquad I=\int\limits_0^1 f(x)\,dx$$

  • a too high order approximation can be less optimal if the function being integrated is not of high order
$$ f(x)=x^2\,,\qquad I=\int\limits_{-1}^2 f(x)\,dx$$

Comparison table

This table shows the number of function evaluations and the asymptotic behaviour of the integration error (under appropriate smoothness condition) for the methods we have seen in this lecture.

method # function evaluation error
Rectangle $N$ $N^{-1}$
Mid-point $N$ $N^{-2}$
Trapezium $N+1$ $N^{-2}$
Simpson $2N+1$ $N^{-4}$
Quartic $4N+1$ $N^{-6}$
6th order $6N+1$ $N^{-8}$

N.B. This is for single variable integrals, it gets worse as the number of variables increases.

Assignment 2

The second assignment which investigates numerical integration using Simpson's rule will be released Monday 2 pm (deadline of the first).