Numerical Integration

Problem: You have a function f and you want to find \int_a^b f(x) dx

Trapezoid Method

\boxed{ A_i = \frac{f(x_i) + f(x_{i-1})}{2} (x_i - x_{i-1}) } \\~\\ \small\textit{(sum for the answer ($A_1 + A_2 + .. + A_n$))}

Example: Applying Trapezoid Method

Q:

f(x) = x^2 + 3x + 1

x_0 = 1, x_1 = 2

Find \int f(x)dx from 1 to 2.

A:

SUPER IMPORTANT NOTE — The subinterval behavior shown in this example (n = 1 \to 2 \to 4 \to 8 \to ...) is not a property of the Trapezoid method. The Trapezoid method itself works for any n.

it1: n=1

A_0 = \left[ \frac{f(1) + f(2)}{2} \right] (1) = 8

it2: n=2

A_1 = \left[ \frac{f(1) + f(1.5)}{2} \right] (0.5) = 3.1875

A_2 = \left[ \frac{f(1.5) + f(2)}{2} \right] (0.5) = 4.6875

A_1 + A_2 = 7.875

it3: n=4

A_1 = \left[ \frac{f(1) + f(1.25)}{2} \right] (0.25) = 1.4140625

A_2 = \left[ \frac{f(1.25) + f(1.5)}{2} \right] (0.25) = 1.7578125

A_3 = \left[ \frac{f(1.5) + f(1.75)}{2} \right] (0.25) = 2.1328125

A_4 = \left[ \frac{f(1.75) + f(2)}{2} \right] (0.25) = 2.5390625

A_1 + A_2 + A_3 + A_4 = 7.8437

Note — This cannot solve non linear functions, unless you do this \infin times.

Romberg’s Method

\boxed{ \text{Romberg's} \begin{cases} R(0,0) &= \frac{f(a) + f(b)}{2} (b - a) \\~\\ R(i,0) &= \frac{R(i-1,0)}{2} + h \sum_{k=1}^{2^{i-1}} f(a + (2k-1) h_i ) \\~\\ R(i,j) &= R(i,j-1) + \frac{R(i,j-1) - R(i-1, j-1)}{4^j - 1} \end{cases} } \\ \small\textit{where $h_i = \frac{h}{2^i}$}

Algorithm: \begin{array}{l | c c c c l} & \mathbf{j=0} & \mathbf{j=1} & \mathbf{j=2} & \cdots & \mathbf{j=N} \\ \hline \mathbf{i=0} & R(0,0) \\ \mathbf{i=1} & R(1,0) & R(1,1) \\ \mathbf{i=2} & R(2,0) & R(2,1) & R(2,2) \\ \mathbf{i=3} & R(3,0) & R(3,1) & R(3,2) & \ddots \\ \vdots & \vdots & \vdots & \ddots & \ddots \\ \mathbf{i=N} & R(N,0) & R(N,1) & \cdots & \cdots & R(N,N) \end{array}

Example: Reading/Constructing the Table

Column zero (j=0) is constructed by the R(0,0) and R(i,0) formulae, which is just the previous trapezoid method and a recursive version of it (the new recursive function just does A_1 + A_2 + ... + A_n in less steps, it isn’t an extrapolation).

The other columns (j>0) is Richardson’s extrapolation, and function the exact same way. R(i,j) depends on the value to its immediate left and diagonally up-left.

The diagonal elements (R(i,i)) are the most-accurate approximations at each iteration i.

Example: Applying Romberg’s

Q:

f(x) = x^2 + 3x + 1

x_0 = 1, x_1 = 2

Find \int f(x)dx from 1 to 2.

A:

it1: h=1

\begin{aligned} R(0,0)&=\frac{f(1)+f(2)}{2} (1)\\ &= 8 \end{aligned}

it2: h=0.5

\begin{aligned} R(1,0) &= \frac{R(0,0)}{2} + (0.5) \sum_{k=1}^{1} f(1 + 0.5) \\ &= \frac{8}{2} + (0.5)f(1.5) \\ &= 7.875 \end{aligned}

Notice — We got the same number as last time, but didn’t need to plug into the function four times.

\begin{aligned} R(1,1) &= R(1,0) + \frac{R(1,0) - R(0,0)}{3} \\ &= 7.875 + \frac{7.875 - 8}{3} \\ &= 7.8\overline{3} \end{aligned}

We’ve converged to the actual value of the integral, something which was impossible with the trapezoid method.

it3: h=0.25

\begin{aligned} R(2,0) &= \frac{R(1,0)}{2} + (0.25)(f(1 + 0.25) + f(1 + (3)(0.25))) \\ &= \frac{7.875}{2} + (0.25)(f(1.25) + f(1.75)) \\ &= 7.84375 \end{aligned}

While the work increased, it’s way less than the 8 plug-ins we did in trapezoidal.

\begin{aligned} R(2,1) &= 7.84375 + \frac{7.84375 - 7.875}{3} \\ &= 7.8\overline{3} \end{aligned}

\begin{aligned} R(2,2) &= 7.8\overline{3} + \frac{7.8\overline{3} - 7.8\overline{3}}{15} \\ &= 7.8\overline{3} \end{aligned}

As you can see, we’ll just continue to get the same answer over-and-over now.

Professor’s Note — For this algorithm, students always falter with the summation.