Numerical Methods/Numerical Integration

<< Numerical Methods

Often, we need to find the integral of a function that may be difficult to integrate analytically (ie, as a definite integral) or impossible (the function only existing as a table of values).

Some methods of approximating said integral are listed below.

Trapezoidal Rule

Consider some function, possibly unknown, $f(x)$, with known values over the interval [a,b] at n+1 evenly spaced points x_i of spacing $h=\frac{(b-a)}{n}$, $x_0=a$ and $x_n=b$.

Further, denote the function value at the ith mesh point as $f(x_i)$.

Using the notion of integration as "finding the area under the function curve", we can denote the integral over the ith segment of the interval, from $x_{i-1}$ to $x_i$ as:

$\underset{x_{i-1}}{\overset{x_i}{\int }}f(x)dx$ = (1)

Since we may not know the antiderivative of $f(x)$, we must approximate it. Such approximation in the Trapezoidal Rule, unsurprisingly, involves approximating (1) with a trapezoid of width h, left height $f(x_{i-1})$, right height $f(x_i)$. Thus,

(1) $\simeq \frac{1}{2}h(f(x_{i-1})+f(x_i))$ = (2)

(2) gives us an approximation to the area under one interval of the curve, and must be repeated to cover the entire interval.

For the case where n = 2,

$\underset{x_a}{\overset{x_b}{\int }}f(x)dx\simeq \frac{1}{2}h(f(x_0)+f(x_1))+\frac{1}{2}h(f(x_1)+f(x_2))$ = (3)

Collecting like terms on the right hand side of (3) gives us:

$\frac{1}{2}h(f(x_0)+f(x_1)+f(x_1)+f(x_2))$

or

$\frac{1}{2}h(f(x_0)+2f(x_1)+f(x_2))$

Now, substituting in for h and cleaning up,

$\frac{(b-a)}{2\cdot 2}(f(x_0)+2f(x_1)+f(x_2))$

To motivate the general version of the trapezoidal rule, now consider n = 4,

$\underset{x_a}{\overset{x_b}{\int }}f(x)dx\simeq \frac{1}{2}h(f(x_0)+f(x_1))+\frac{1}{2}h(f(x_1)+f(x_2))+\frac{1}{2}h(f(x_2)+f(x_3))+\frac{1}{2}h(f(x_3)+f(x_4))$

Following a similar process as for the case when n=2, we obtain

$\frac{(b-a)}{2\cdot 4}(f(x_0)+2(f(x_1)+f(x_2)+f(x_3))+f(x_4))$

Proceeding to the general case where n = N,

$\underset{x_a}{\overset{x_b}{\int }}f(x)dx\simeq \frac{(b-a)}{2\cdot n}(f(x_0)+2(\sum _{k=1}^{N}f(x_k))+f(x_n))$

This is an example of what the trapezoidal rule would represent graphicly, here $y=-x^2+5$.

Example

Approximate $\underset{0}{\overset{1}{\int }}{x}^{3}dx$ to within 5%.

First, since the function can be exactly integrated, let us do so, to provide a check on our answer.

$\underset{0}{\overset{1}{\int }}{x}^{3}dx={\left[\frac{x^4}{4}\right]}_0^1=\frac{1}{4}=0.25$ = (4)

We will start with an interval size of 1, only considering the end points.

$f(0)=0$

$f(1)=1$

(4) $\simeq \frac{(1-0)}{(2\cdot 1)}(f(0)+f(1))=\frac{1}{2.1}(0+1)=\frac{1}{2}=0.5$

Relative error = $\left|\frac{(0.5-0.25)}{0.25}\right|=1$

Hmm, a little high for our purposes. So, we halve the interval size to 0.5 and add to the list

$f(0.5)=0.125$

(4) $\simeq \frac{(1-0)}{(2\cdot 2)}(f(0)+2f(0.5)+f(1))=\frac{1}{2\cdot 2}(0+2(0.125)+1)=\frac{1.25}{4}=0.3125$

Relative error = $\left|\frac{(0.3125-0.25)}{0.25}\right|=0.25$

Still above 0.01, but vastly improved from the initial step. We continue in the same fashion, calculating $f(0.25)$ and $f(0.75)$, rounding off to four decimal places.

$f(0.25)=0.0156$

$f(0.75)=0.4219$

(4) $\simeq \frac{(1-0)}{(2\cdot 4)}(0+2(0.0156+0.125+0.4219)+1)=\frac{1}{8}(2.2150)=0.2656$

Relative error = $\left|\frac{(0.2656-0.25)}{0.25}\right|=0.0624$

We are well on our way. Continuing, with interval size 0.125 and rounding as before,

$f(0.125)=0.0020$

$f(0.375)=0.0527$

$f(0.625)=0.2441$

$f(0.875)=0.6699$

(4) $\simeq \frac{(1-0)}{(2\cdot 8)}(0+2(0.0020+0.0156+0.0527+0.0125+0.2441+0.4219+0.6699)+1)=\frac{1}{16}(4.0624)=0.2539$

Relative error = $\left|\frac{(0.2539-0.25)}{0.25}\right|=0.0156$

Since our relative error is less than 5%, we stop.

Error Analysis

Let y=f(x) be continuous,well-behaved and have continuous derivatives in [x₀,x_n]. We expand y in a Taylor series about x=x₀,thus-
$\underset{x_0}{\overset{x_1}{\int }}ydx=\underset{x_0}{\overset{x_1}{\int }}[y_0+(x-x_0)y{\text{'}}_0+(x-x_0{)}^2{y}^{\prime }{\text{'}}_0/2!+......]dx$

Simpson's Rule

Consider some function $y=f(x)$ possibily unknown with known values over the interval [a,b] at n+1 evently spaced points then it defined as

$\underset{x_0}{\overset{x_n}{\int }}f(x)dx\simeq \frac{1}{3}h\{f(x_0)+f(x_n)+2(f(x_2)+f(x_4)+...+f(x_{n-2}))+4(f(x_1)+f(x_3)+...+f(x_{n-1})\left)\right\}$

where $h=\frac{(b-a)}{n}$ and $x_0=a$ and $x_n=b$.

Example

Evaluate $\underset{0}{\overset{1.2}{\int }}x{\left(8-x^3\right)}^{\frac{1}{2}}dx$ by taking $n=6$ ($n$ must be even)

Solution: Here $f(x)=x{\left(8-x^3\right)}^{\frac{1}{2}}$

Since $a=0$ & $b=1.2$ so $h=\frac{b-a}{n}=\frac{1.2-0}{6}=0.2$

Now when $a=x_0=0$ then $f(x_0)=0$

And since $x_n=x_{n-1}+h$, therefore for $x_1=0.2$ , $x_2=0.4$ , $x_3=0.6$ , $x_4=0.8$ , $x_5=1$ , $x_6=b=1.2$ the corresponding values are $f(x_1)=0.7784$ , $f(x_2)=1.58721$ , $f(x_3)=1.6740$ , $f(x_4)=2.1891$ , $f(x_5)=2.6458$ , $f(x_6)=3.0053$

Incomplete ... Completed soon

References

Eric W. Weisstein. "Trapezoidal Rule." From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/TrapezoidalRule.html

---

Main Page - Mathematics bookshelf - Numerical Methods