# Numerical Integration
## CH EN 2450 - Numerical Methods
**Prof. Tony Saad (<a>www.tsaad.net</a>) <br/>Department of Chemical Engineering <br/>University of Utah**
<hr/>

The purpose of numerical integration is to compute the area under a curve or a specified set of data points. The two general applications of numerical integration comprise the integration of (a) complex continuous functions, and (b) discrete data. Many times, the function to be integrated does not admit an antiderivative, such as $f(x) = e^{-x^2}$.

## Newton-Cotes Integration
The general idea behind Newton-Cotes Integration is to approximate the function using a polynomial and then integrating the polynomial.

### Left, Right, and Midpoint Rules
The left, right, and midpoint rules are the simplest of all numerical integration methods. They are based on approximating the function as a constant
\begin{equation}
\int_a^b f(x) \text{d}x \approx \int_a^b C \text{d}x = C\times(b-a)
\end{equation}
The constant can be chosen in several ways, the most common are
1. Left point rule: choose $C = f(a)$
2. Right point rule: choose $C = f(b)$
3. Midpoint rule: choose $C = f(\frac{a+b}{2})$

The three options are shown on the plots below.


```python
import numpy as np
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
import matplotlib.pyplot as plt
```


```python
def myf(x):
    return 1.0/np.sqrt(np.pi) * np.exp(-x**2)
```


```python
xalot = np.linspace(0,2,200)

npts = 5

x = np.linspace(0,2, npts)
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, sharey=True)
fig.set_size_inches((12,3))
a = 0.5
b = 1.0

ax = ax1
ax.plot(xalot, myf(xalot))
ax.set_title('Leftpoint Rule')
# ax1.fill_between([a,b],[myf(a),myf(a)],  edgecolor='k', facecolor='lightgray',linewidth=1)
# plot left Riemann Sum
for i in range(0, len(x)-1):
    a = x[i]
    b = x[i+1]
    ax.fill_between([a,b],[myf(a),myf(a)],  edgecolor='k', facecolor='lightgray',linewidth=1)

ax = ax2
ax.plot(xalot, myf(xalot))
# ax3.fill_between([a,b],[myf( (a+b)/2) ,myf( (a+b)/2)],  edgecolor='k', facecolor='lightgray',linewidth=1)
# plot midpoint rule
ax.set_title('Midpoint Rule')
for i in range(0, len(x)-1):
    a = x[i]
    b = x[i+1]
    h = myf((a+b)/2)
    ax.fill_between([a,b],[h,h],edgecolor='k', facecolor='lightgray',linewidth=1)

ax = ax3
ax.plot(xalot, myf(xalot))
# ax2.fill_between([a,b],[myf(b),myf(b)],  edgecolor='k', facecolor='lightgray',linewidth=1)
# plot right Riemann Sum
ax.set_title('Rightpoint Rule')
for i in range(0, len(x)-1):
    a = x[i]
    b = x[i+1]
    ax.fill_between([a,b],[myf(b),myf(b)],  edgecolor='k', facecolor='lightgray',linewidth=1)

plt.ylim(bottom=0.0)
```




    (0.0, 0.5923990627251441)




    

    


Here we define all three rules


```python
def leftpoint(f, a, b, npts):
    '''
    f: Any Python function
    a: Lower integral bound
    b: Upper integral bound
    npts: Number of quadrature points
    
    Returns the integral of f(x) based on the leftpoint rule
    '''    
    x = np.linspace(a,b,npts)
    sum = 0.0
    for i in range(0,len(x)-1):
        a = x[i]
        b = x[i+1]
        sum += (b-a)*(f(a) )
    return sum

def rightpoint(f, a, b, npts):
    '''
    f: Any Python function
    a: Lower integral bound
    b: Upper integral bound
    npts: Number of quadrature points
    
    Returns the integral of f(x) based on the rightpoint rule
    '''
    x = np.linspace(a,b,npts)
    sum = 0.0
    for i in range(0,len(x)-1):
        a = x[i]
        b = x[i+1]
        sum += (b-a)*(f(b) )
    return sum


def midpoint(f, a, b, npts):
    '''
    f: Any Python function
    a: Lower integral bound
    b: Upper integral bound
    npts: Number of quadrature points
    
    Returns the integral of f(x) based on the midpoint rule
    '''
    x = np.linspace(a,b,npts)
    sum = 0.0
    for i in range(0,len(x)-1):
        a = x[i]
        b = x[i+1]
        sum += (b-a)*(f( (a+b)/2 ) )
    return sum
```

#### Example

Use the left, right, and midpoints rules to compute the integral 
\begin{equation}
\int_0^3 e^{-x^2} \text{d} x
\end{equation}
Compute the absolute relative approximate error for a set of integration points ranging from 2 to 200. Plot the error in all cases.


```python
numpts = np.arange(2,50)
a = 0.0
b = 3.0

lpvals = []
rpvals = []
mpvals = []
for npts in numpts:
    v = leftpoint(myf,a,b,npts)     
    lpvals.append(v)
    v = rightpoint(myf,a,b,npts)
    rpvals.append(v)
    v = midpoint(myf,a,b,npts)
    mpvals.append(v)
    
plt.plot(numpts,lpvals, numpts,rpvals, numpts, mpvals)
plt.grid()
```


    

    


Now let's look at how the absolute approximate relative error behaves


```python
lpvals = np.array(lpvals)
e1 =abs( (lpvals[1:] - lpvals[:-1])/lpvals[1:] )
plt.loglog(numpts[1:],e1, label="Leftpoint rule")

rpvals = np.array(rpvals)
e2 =abs( (rpvals[1:] - rpvals[:-1])/rpvals[1:])
plt.loglog(numpts[1:],e2, label="Rightpoint rule" )

mpvals = np.array(mpvals)
e3 =abs( (mpvals[1:] - mpvals[:-1])/mpvals[1:])
plt.loglog(numpts[1:],e3, label="Midpoint rule" )

plt.grid()
plt.legend()
```




    <matplotlib.legend.Legend at 0x11b2990b8>




    

    


### Trapezoid Rule

The trapezoidal rule fits a straight line between each two points. It is given by the formula
\begin{equation}
\int_a^b f(x) \text{d} x \approx \frac{b-a}{2}[f(a) + f(b)]
\end{equation}


```python
def myf(x):
    return 1.0/np.sqrt(np.pi) * np.exp(-x**2)

xalot = np.linspace(0,2,200)

npts = 5

x = np.linspace(0,2, npts)

plt.plot(xalot, myf(xalot))
plt.title('Trapezoidal Rule')
# ax1.fill_between([a,b],[myf(a),myf(a)],  edgecolor='k', facecolor='lightgray',linewidth=1)
# plot left Riemann Sum
for i in range(0, len(x)-1):
    a = x[i]
    b = x[i+1]
    plt.fill_between([a,b],[myf(a),myf(b)],  edgecolor='k', facecolor='lightgray',linewidth=1)

plt.ylim(bottom=0.0)
```




    (0.0, 0.5923990627251441)




    

    


Implementation of the Trapezoidal rule is straightforward as shown below. In this case, we will implement the trapezoidal rule for both continuous functions and for discrete datasets.


```python
def traprule(x,y):
    '''
    Integrates a discrete set of data using the Trapezoidal rule.
    
    x: Values of the independent variable
    y: Values of the dependent variable
    
    Returns the integral of f(x) based on the midpoint rule
    '''

    sum = 0.0
    for i in range(0,len(x)-1):
        sum += 0.5*(x[i+1]-x[i])*(y[i+1] + y[i])
    return sum


def traprulef(f, a, b, npts):
    '''
    Computes the integral of f(x) using the Trapezoidal rule. 
    
    f: Any Python function
    a: Lower integral bound
    b: Upper integral bound
    npts: Number of quadrature points
    
    Returns the integral of f(x) based on the midpoint rule
    '''    
    x = np.linspace(a,b,npts)
    sum = 0.0
    for i in range(0,len(x)-1):
        a = x[i]
        b = x[i+1]
        sum += 0.5*(b-a)*(f(a) + f(b))
    return sum
```

#### Example
Use the Trapezoidal rule to compute the integral 
\begin{equation}
\int_0^3 e^{-x^2} \text{d} x
\end{equation}
Compute the absolute relative approximate error for a set of integration points ranging from 2 to 200. Plot the error.


```python
numpts = np.arange(2,30)
a = 0.0
b = 3.0

trapvals = []
for npts in numpts:
    v = traprulef(myf,a,b,npts)     
    trapvals.append(v)
    
plt.plot(numpts,trapvals)
plt.grid()
```


    

    



```python
trapvals = np.array(trapvals)
e1 =abs( (trapvals[1:] - trapvals[:-1])/trapvals[1:] )
plt.loglog(numpts[1:],e1, label="Trapezoidal rule")
plt.grid()
```


    

    



```python
from IPython.core.display import HTML
def css_styling():
    styles = open("../../styles/custom.css", "r").read()
    return HTML(styles)
css_styling()
```




CSS style adapted from https://github.com/barbagroup/CFDPython. Copyright (c) Barba group
<link href='http://fonts.googleapis.com/css?family=Merriweather' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Bitter' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Oxygen' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Lora' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>

<style>
    @font-face {
        font-family: "Computer Modern";
        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');
    }

    /*div.cell{
        width:800px;
        margin-left:16% !important;
        margin-right:auto;
    } */

    /* set the font size in tables */
    tr, td, th{
        font-size:110%;
    }

    /* spec for headers */
    h1 {
        font-family: 'Bitter', serif;
    }

    h2 {
        font-family: 'Fenix', serif;
    }

    h3{
        font-family: 'Fenix', serif;
        margin-top:12px;
        margin-bottom: 3px;
    }

    h4{
        font-family: 'Fenix', serif;
    }

    h5 {
        font-family: 'Alegreya Sans', sans-serif;
    }

    div.text_cell_render{
        font-family: 'Merriweather','Alegreya Sans','Lora', 'Oxygen', "Helvetica Neue", Arial, Helvetica, Geneva, sans-serif;
        line-height: 160%;
        font-size: 130%;
    }

    .CodeMirror{
        font-family: "Source Code Pro";
        font-size: 100%;
    }

    .text_cell_render h1 {
        font-weight: 200;
        font-size: 32pt;
        line-height: 120%;
        color:#CD2305;
        margin-bottom: 0.5em;
        margin-top: 0.5em;
        display: block;
    }

    .text_cell_render h2 {
        font-size: 26pt;
        text-align: center;
    }

    .text_cell_render h3 {
        font-size: 20pt;
    }

    .text_cell_render h4 {
        font-size: 18pt;
    }

    .text_cell_render h5 {
        font-weight: 300;
        font-size: 16pt;
        color: #CD2305;
        font-style: italic;
        margin-bottom: .5em;
        margin-top: 0.5em;
        display: block;
    }

    .warning{
        color: rgb( 240, 20, 20 )
    }  

/*  div#notebook {background-color: #1e1e1e; border-top: none;}
    div#notebook-container {background-color: rgb(180, 180, 180);}
 */

</style>







```python

```

f622b9266a9bb5e2524ce3fa8341e6c2c3042293
ipynb
Jupyter Notebook
topics/numerical-integration/numerical integration.ipynb
jomorodi/NumericalMethods
e040693001941079b2e0acc12e0c3ee5c917671c
MIT
2019-03-27T05:22:34.000Z
2021-01-27T10:49:13.000Z
topics/numerical-integration/numerical integration.ipynb
jomorodi/NumericalMethods
e040693001941079b2e0acc12e0c3ee5c917671c
MIT
topics/numerical-integration/numerical integration.ipynb
jomorodi/NumericalMethods
e040693001941079b2e0acc12e0c3ee5c917671c
MIT
2019-12-29T23:31:56.000Z
2021-12-28T19:04:10.000Z
Qwen/Qwen-72B
1. YES
2. YES