📂Numerical Analysis

Numerical Integration

Definition ¹

Let’s assume that $f : [a,b] \to \mathbb{R}$ is integrable over $[a,b]$, and $[a,b]$ is divided into nodes like $a = x_{0} < \cdots < x_{n} = b$.

1. The integral operator $I$ is defined as $\displaystyle I(f) := \int_{a}^{b} f(x) dx$.
2. The integral operator $I_{n}$ is defined as $\displaystyle I_{n} (f) := \sum_{k=1}^{n} \int_{x_{k-1}}^{x_{k}} f(x) dx$.
3. The error $E_{n}$ is defined as $E_{n} (f) := I (f) - I_{n} ( f )$.
4. $\displaystyle \lim_{n \to \infty} {{\tilde{E}_{n} (f) } \over { E_{n} (f) }} = 1$, satisfying $\tilde{E}_{n}$, is referred to as the asymptotic error for $E_{n}$.

Explanation

No one would dispute that calculating definite integrals is one of the goals of numerical analysis. In fact, the idea behind the method of exhaustion suggests that the more we divide the space, the better we can approximate the value of the integral. Thus, when it comes to numerical integration of interest in numerical analysis, it is about how much error there is, how to integrate when not in a closed interval, etc.

Example

Trapezoidal Rule

As an example of numerical integration, let’s take a look at the trapezoidal rule for $f \in C^2 [a,b]$. The Trapezoidal Rule approximates the value for $I(f)$ by linearly interpolating $f$ over the interval $[a,b]$ and taking the integral of that function as the approximation for $I(f)$.

$$I_{1}^{1} (f) := \left( {{ b - a } \over { 2 }} \right) [ f(a) + f(b) ]$$ In this case, the actual error between $I(f)$ and $I_{1}^{1} (f)$ $E_{1}^{1} (f)$ can be computed as follows for some $\xi \in [a,b]$.

Polynomial Interpolation:

• [4] Error with the actual function: For a differentiable $f : \mathbb{R} \to \mathbb{R}$ up to $(n+1)$ times and for some $\xi \in \mathscr{H} \left\{ x_{0} , \cdots , x_{n} \right\}$, the polynomial interpolation of $f$ $p_{n}$ satisfies $\displaystyle f(t) - p_{n} (t) = {{ (t - x_{0}) \cdots (t - x_{n}) } \over { (n+1)! }} f^{(n+1)} ( \xi )$ for some $t \in \mathbb{R}$.

$$\begin{align*} \displaystyle E_{1}^{1} (f) :=& I(f) - I_{1}^{1} (f) \\ =& \int_{a}^{b} \left[ f(x) - {{ f(b) ( x - a ) - f(a) (x - b) } \over { b - a }} \right] dx \\ =& \int_{a}^{b} \left[ f(x) - p_{1} (x) \right] dx \\ =& {{1} \over {2}} f '' ( \xi ) \int_{a}^{b} (x-a) (x-b) dx \\ =& \left[ {{1} \over {2}} f '' ( \xi ) \right] \left[ - {{1} \over {6}} (b-a)^{3} \right] \\ =& - {{1} \over {12}} (b-a)^{3} f '' ( \xi ) \end{align*}$$

Composite Trapezoidal Rule

Now, considering dividing $[a,b]$ into $n$ intervals and approximating the integral at each $[x_{k-1} , x_{k} ]$ using the trapezoidal rule, the method is called the (Composite) Trapezoidal Rule. $$I_{n}^{1} (f) := \sum_{k=1}^{n} {{ h } \over { 2 }} \left( f (x_{k-1}) + f( x_{k} ) \right)$$ Then, the actual error between $I(f)$ and $I_{n}^{1} (f)$ can be computed as follows for some $\xi_{k} \in [x_{k-1}, x_{k} ]$. $$\begin{align*} \displaystyle E_{n}^{1} (f) =& I (f) - I_{n}^{1} (f) \\ =& \sum_{k=1}^{n} \left( - {{ h^3 } \over { 12 }} f '' ( \xi_{k} ) \right) \\ =& - {{ h^3 n } \over { 12 }} \left( {{ 1 } \over { n }} \sum_{k=1}^{n} f '' ( \xi_{k} ) \right) \end{align*}$$ Here, since $\displaystyle \left( {{ 1 } \over { n }} \sum_{k=1}^{n} f '' ( \xi_{k} ) \right)$ is the average of $f’’ ( \xi_{k} )$, due to the Intermediate Value Theorem, there must exist some $\xi \in [a,b]$ that satisfies $\displaystyle \min_{x \in [a,b]} f ''(x) \le f ''( \xi ) \le \max_{x \in [a,b]} f ''(x)$. Thus, $$E_{n}^{1} (f) = - {{ (b-a) h^2 } \over {12}} f '' (\xi )$$ Seeing the error formulated like this lets us understand that the wider the interval $[a,b]$, the larger the error becomes, and the smaller $h$ is, the more accurate it gets. However, due to $f’’ ( \xi )$, we can’t know how much exactly it offsets. To overcome this, the asymptotic error is calculated as follows. $$\begin{align*} \lim_{n \to \infty} {{ E_{n}^{1} (f) } \over { h^2 }} =& \lim_{n \to \infty} {{1} \over {h^2}} \sum_{k=1}^{n} \left( - {{ h^3 } \over { 12 }} f '' ( \xi_{k} ) \right) \\ =& - {{ 1 } \over { 12 }} \lim_{ n \to \infty} \sum_{k=1}^{n} h f '' ( \xi_{k} ) \\ =& - {{ 1 } \over { 12 }} \int_{a}^{b} f ''(x) dx \\ =& - {{ 1 } \over { 12 }} [ f '(b) - f '(a) ] \end{align*}$$ If we only know the error $E_{n}^{1} (f)$, according to the Theorem, we know that $f’’ ( \xi )$ exists but couldn’t specify its value. But now, $$\tilde{E}_{n}^{1} (f) := \lim_{n \to \infty} E_{n}^{1} (f)$$ for sufficiently large $n$, it can be understood that $$E_{n}^{1} (f) \approx - {{ h^2 } \over { 12 }} [ f '(b) - f '(a) ] = \tilde{E}_{n}^{1} (f)$$ which makes it relatively simple and accurate to ascertain the error. Moreover, there’s no reason to stop at just knowing the extent of the error. If we can predict the error before even calculating, we can also correct for the error from the beginning. $$I(f) - I_{n}^{1} (f) \approx \tilde{E}_{n}^{1} (f) \implies I(f) \approx I_{n}^{1} (f) + \tilde{E}_{n}^{1} (f)$$

Corrected Trapezoidal Rule

Let’s define the operator for calculating $I(f)$ as follows. $$CT_{n} (f) := h \left[ {{1} \over {2}} f (x_{0} ) + f (x_{1} ) + \cdots + f (x_{n-1} ) + {{1} \over {2}} f (x_{n} ) \right] - {{ h^2 } \over { 12 }} [ f '(b) - f '(a) ]$$

This method is referred to as the (Corrected) Trapezoidal Rule. [ NOTE: If you’ve followed the improvement process of the trapezoidal rule, you’d realize that there’s no need to distinguish between them explicitly. Simply calling it the trapezoidal rule would convey the meaning. ]

Implementation

The results of calculating $\displaystyle \int_{0}^{\pi} e^{x} \cos (x) dx = -{{e^{\pi} + 1} \over {2}} \approx - 12.0703$ using $I_{n}^{1}$ and $CT_{n}$ are as follows.

The first column represents $n$, the second $\left| I(f) - I_{n}^{1} (f) \right|$, and the third $\left| I(f) - CT (f) \right|$. It can be observed that the error decreases by a factor of $\displaystyle {{1} \over {4}}$ when the number of nodes is doubled in the composite trapezoidal rule, and by $\displaystyle {{1} \over {16}}$ in the corrected trapezoidal rule. It’s evident that the corrected trapezoidal rule performs significantly better, especially in terms of speed.

Below is an example code written in Python.

import numpy as np
def f(x) :
    return np.exp(x)*np.cos(x)
 
def d(f,x,tol=10**(-8)) :
    h=1
    d1=1
    d2=0
 
    while(abs(d1-d2)>tol) :
        d1=d2
        d2=(f(x+h)-f(x))/h
        h=h/2
 
    return(d2)
 
def I(f,a,b,n) :
    h=(b-a)/n
    x=np.linspace(a,b,n+1)
    x=f(x)
    x[0]=x[0]/2
    x[-1]=x[-1]/2
    return(h*sum(x))
 
def CT(f,a,b,n) :
    h=(b-a)/n
    x=np.linspace(a,b,n+1)
    x=f(x)
    x[0]=x[0]/2
    x[-1]=x[-1]/2
    return(h*sum(x)-(h**2)*(d(f,b)-d(f,a))/12)
 
 
TV=-(np.exp(np.pi)+1)/2
print("True Value : %2.4f" % TV)
 
print("|   n |  I(f) Err | CT(f) Err |")
print("="*30)
for n in range(1,10) :
    print("| %3d |  %8.2e |  %8.2e |" % (2**n,abs(TV-I(f,0,np.pi,2**n)),abs(TV-CT(f,0,np.pi,2**n))))