Numerical Integration:

The problem of numerical integration is that of determining an approximate value of the integral

$$I=\int\limits_{a}^{b}f(x)dx$$

If $p(x)$ is the interpolating polynomial of degree n which approximate $f(x)$ on $[a,b]$ then for the error

$$E=\int\limits_{a}^{b}f(x)dx-\int\limits_{a}^{b}p(x)dx$$

we have the following theorem:
Theorem: Let $x_i \,(i=0,1,・,n)$ be $n+1$ points on the interval $[a,b]$. Let $f(x_i)$ be the corresponding values of $f$ at $x=x_i$. Let $p(x)$ be the polynomial of degree n which interpolates at the n+1 points $(x_i,f_i)\,(i=0,1,...n).$ Then if

$$\psi(x)=(x-x_0)(x-x_1)...(x-x_n)$$

does not change sign on the interval $[a,b]$, the error of numerical integration is given by

$$\int\limits_{a}^{b}f(x)dx-\int\limits_{a}^{b}p(x)dx=\frac{f^{(n+1)}(\eta)}{(n+1)!}\int\limits_{a}^{b}\psi(x)dx$$

For some point $\eta$ in $[a,b]$.
Proof: We have from the error formula for the interpolation polynomial

$$f(x)-p(x)=\frac{\psi(x)f^{(n+1)}(\xi)}{(n+1)!}$$

Integration, we get

$$\int\limits_{a}^{b}f(x)dx-\int\limits_{a}^{b}p(x)dx=\int\limits_{a}^b\frac{\psi(x)f^{(n+1)}(\xi)dx}{(n+1)!}$$

If $\psi(x)$ does not change sign on the interval $[a,b]$, then we can apply the second mean value theorem of integral calculus, we obtain the desired result. Let us assume that the function $f(x)$ can be computed at a set of n+1 equally spaced points $x_0, x_1,...,x_n$. Then $f(x)$ can be approximated by the Newton forward difference interpolating polynomial

$$f(x)\simeq p(s)=f_0+S\Delta f_0+\left(% \begin{array}{c} S... ...gin{array}{c} S \\ n \\ \end{array}% \right)\Delta^nf_n$$

where $S=(x-x_0)/h$. Thus we have

$$\int\limits_{a}^{b}f(x)dx\simeq h\int\limits_0^1P(S)dS$$

Let us take n=0, we get

$$\int\limits_{a}^{b}f(x)dx\simeq h\int\limits_0^1f_0dS=hf_0$$

This is known as rectangular rule and we denote it as

$$R_1=h f_0$$

The error of this approximation is

$$E_1^R=f'(\eta_0)\int\limits_{x_0}^{x_1}(x-x_0)dx=h^2\frac{f'(\eta_0)}{2}\qquad x_0 < \eta_0 x_1$$

To find the integral of f(x) over an extended interval $[a,b]$ we subdivide $[a,b]$ into N equal subdivision, setting $x_0=a, x_1=h+h,...x_N=h+Nh=b$ Now

$$\int\limits_{a}^{b}f(x)dx=\int\limits_{a}^{x_1}f(x)dx+\int\limits_{x_1}^{x_2}f(x)dx+...+\int\limits_{x_{N-1}}^{b}f(x)dx$$

Applying the above formula to each integral of yields the rectangular rule for the integral of $f(x)$ over an interval $[a,b]$:

$$R_N=h(f_0+f_1+...+f_{N-1})$$

and the error is given by

$$E^R_N=\frac{ h^2}{2}[f'(\eta_0)+f'(\eta_1)+...+f'(\eta_{N-1})]$$

and if $f('x)$ is continuous over $[a,b]$

$$E^R_N=\frac{(b-a)}{2}hf'(\eta)\qquad a < \eta < b$$

A more accruable formula can be obtained by taking n=1, and we find that

$$\int\limits_{a}^{b}f(x)dx=hh\int\limits_{0}^1(f_0+S\Delta f_0)dS=\frac{h}{2}(f_0+f_1)$$

which is known as Trapezoidal rule denoted by

$$T_1=\frac{h}{2}(f_0+f_1)$$

with the error given by

$$E^T_1=f''\frac{\eta_0}{2}\int\limits_{a}^{b}(x-x_0)(x-x_1)dx=-\frac{h^3f''(\eta_0)}{12}$$

To obtain the integral over the interval $[a,b]$, we subdivide $[a,b]$ into N equal parts and use the above formula to each integral, this yields

$$T_N=h(\frac{f_0}{2}+f_1+f_2+...+f_{N-1}+\frac{f_N}{2})$$

The error of this formula is given by

$$E^T_N=-\frac{(b-a)}{12}h^2f''(\eta) \qquad 0 < \eta < b$$

Simpson's Rule
We integrate $f(x)$ over the double interval $[x_0, x_2]$ of width $2h$ and get

$$\int\limits_{x_0}^{x_2}f(x)dx=h\int\limits_{0}^{2}P(S)dS=\int... ...ray}{c} S \\ 4 \\ \end{array}% \right)\Delta ^4f_0+..dS$$

By direct integration, we find that

$$\int\limits_{0}^{2}dS=2,\, \int\limits_{0}^{2}SdS=2,\,\int\li... ...array}{c} S \\ 4 \\ \end{array}% \right)dS=\frac{1}{90}$$

If we retain difference through the third order, we obtain an approximation

$$S_2=h(2f_0+2\Delta f_0+\frac{1}{2}\Delta ^2 f_0+0.\Delta^3f_0)$$

$$=\frac{h}{3}(f_2)+4f_1+f_0$$

This formula is called Simpson's rule.
The error of this formula is given by

$$E^S_2=-\frac{h^5}{90}f^{IV}(\eta)\qquad x_0 < n < x_2$$

To extend Simpson's rule in integration over an interval $[a,b]$, we now divide $[a,b]$ into an even number 2N of sub intervals of width h so that

$$x_0=a\qquad x_{2N}=b \qquad (b-a)=2Nh$$

and

$$x_i-x_{i-1}=h\qquad i=1,2,...2N$$

Using Simpson's rule over the interval $[x_i,x_{i+2}]\,(i=0,2,...,2N-2)$,we have

$$\int\limits_{x_i}^{x_{i+2}}f(x)dx=\frac{h}{2}(f_{i+2}+4f_{i+1}+f_i)-\frac{h^5}{90}f^{IV}(\eta_i)\qquad x_i < \eta_i < x_{i+2}$$

If we now sum over the N subgroups of the intervals each, we obtain

$$\int\limits_{a}^{b}f(x)dx=\frac{h}{3}(f_0+4f_1+2f_2+...+4f_{2N-1}+f_{2N})-\frac{N}{90}h^5f^{IV}{\eta}$$

and using Simpson's rule for integration over an interval $[a,b]$ which has been subdivided into 2N subintervals of length h is

$$S_{2N}=\frac{h}{3}(f_0+4f_1+2f_2+4f_3+...+4f_{2N-1}+f_{2N})$$

and since $N=\frac{b-a}{2h}$, the error term is

$$E^S_{2N}=-\frac{(b-a)}{180}h^4f^{IV}(\eta) \qquad a < \eta < b.$$