<a href="https://colab.research.google.com/github/cohmathonc/biosci670/blob/master/IntroductionComputationalMethods/04_IntroCompMethods_NumericalIntegration.ipynb" target="_parent"></a>


```python
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
```

## Numerical Integration

Numerical integration methods seek an approximate solution to a definite integral
$$\int_a^b f(x)\,  dx \, .$$


```python
# define function
def f(x):
    return 1/(1+np.power(x,2))

# domain bounds
a = 0
b = 3

# plot
fig = plt.figure(figsize=plt.figaspect(0.5))
ax = fig.add_subplot(111) 

x = np.linspace(a-0.1*(b-a), b+0.1*(b-a), 100)
y = f(x)
ax.plot(x, y,'-', label='$f(x)$')

X = np.linspace(a,b,100)
Y = f(X)
ax.fill_between(X,Y, color='r', alpha=0.2, label="$\int_a^b f(x)\,  dx$")

ax.legend(prop={'size':15})
plt.show()

```

In above example, we have chosen 
$$f(x)=\frac{1}{1+x^2} \; .$$
The integral of this function is known:
$$\int \frac{1}{1+x^2} \, dx= \tan^{-1}(x) + C =\arctan(x) + C$$
So we can solve the integral analytically:
$$\int_a^b \frac{1}{1+x^2} \, dx = \arctan(b) - \arctan(a)$$



```python
# plot arctan(x)
fig = plt.figure(figsize=plt.figaspect(0.5))
ax = fig.add_subplot(111) 

y = np.arctan(x)
ax.plot(x, y,'-', label='arctan(x)')
ax.axvline(x=a, color='green', linestyle=':')
ax.axvline(x=b, color='red', linestyle=':')
ax.axhline(y=np.arctan(a), color='green', linestyle='--', label='arctan(x=a)')
ax.axhline(y=np.arctan(b), color='red', linestyle='--',  label='arctan(x=b)')

ax.legend(prop={'size':15})
plt.show()
```

However, not all functions can be integrated analytically, and even if a closed formulation exists, it may be easier to compute a numerical approximation rather than the antiderivative. Also, the integrand $f(x)$ may only be known at certain points.

### Riemann Sums


The [Riemann Integral](https://en.wikipedia.org/wiki/Riemann_integral) is a natural starting point for considering approximation methods for integrals.

Let's assume that $f(x)$ is a bounded function defined on the interval $[a, b]$ which is divided into $N$ subintervals of length $\Delta x = \left| a-b\right|/N$ defining a partition $\{x_0, x_1, \ldots x_{N}\}$ where $a=x_0 < x_1 < x_2 < \ldots <x_{N-1} < x_{N} =b$. 

In each of the $N$ subintervals $[x_i, x_{i+1}]\; i=0, 1, \ldots, N-1$, the function $f(x)$ will have a maximum value $f_{\text{max, i}}$ and a minimum value $f_{\text{min, i}}$.
Defining $\Delta_i=x_{i+1}-x_i$, we can compute an upper and lower limit for the "area under the curve" $\mathcal{I}$ by summing the contribution of all subintervals:

\begin{align}
\text{lower bound:}& \qquad \mathcal{s}_N  =  \sum_{i=0}^{N-1} f_{\text{min, i}}\, \Delta_i \leq \mathcal{I} \tag{1a}\\
\text{upper bound:}& \qquad \mathcal{S}_N =  \sum_{i=0}^{N-1} f_{\text{max, i}}\, \Delta_i \geq \mathcal{I}\tag{1b}\\
\end{align}

These are called the *lower* and *upper* [Riemann Sum](https://en.wikipedia.org/wiki/Riemann_sum#Trapezoidal_rule), respectively.
As $N$ approaches infinity, i.e. $\Delta_i\rightarrow 0$, $\mathcal{s}_N$ and $\mathcal{S}_N$ converge to the exact value $\mathcal{I}$.



#### Rectangular Method: Left & Right Sums

Definitions (1a), (1b) can be used to approximate an integral numerically.
Instead of using the minimum or maximum of the function in each subinterval, we can use the value of the function on either side of each subinterval: 

\begin{align}
\text{left sum:}& \qquad \mathcal{S}_{L, N}  =  \sum_{i=0}^{N-1} f(x_i)\, \Delta_i \approx \mathcal{I} \tag{2a}\\
\text{right sum:}& \qquad \mathcal{S}_{R, N} =  \sum_{i=0}^{N-1} f(x_{i+1})\, \Delta_i \approx \mathcal{I}\tag{2b}\\
\end{align}

These methods are examples of the **rectangular rule** because the area corresponding to each subinterval is approximated by a rectangle with dimensions length of the subinterval $\Delta_i$ and the function value at one of the endpoints of the subinterval.

The **left rule** will overestimate the actual value of the integral $\mathcal{I}$ in a region where f(x) decreases monotonically, and will overestimate its value where f(x) is increase monotonically. The opposite holds for the **right rule**.


#### Rectangular Method: Midpoint Rule

Another variation of the rectangular method is the **midpoint rule** , where the function value at the midpoint (instead of one of the endpoints) is used to compute the  rectangular area corresponding to each subinterval:


\begin{align}
\text{midpoint sum:}& \qquad \mathcal{S}_{M, N}  =  \sum_{i=0}^{N-1} f\left(\frac{x_i + x_{i+1}}{2}\right)\, \Delta_i \approx \mathcal{I} \tag{3}\\
\end{align}

As the comparison below shows, the midpoint rule provides a more accurate approximation of the integral than the rectangular method.

#### Comparison: Left, Right, Midpoint Rule


```python
# define function
def f(x):
    return 1/(1+np.power(x,2))

# domain bounds
a = 0
b = 3

# actual function, evaluated on finer grid to be represented smoothly
X = np.linspace(a-0.1*(b-a), b+0.1*(b-a), 100)
Y = f(X)

# partition for numeric approximation
N = 10
x = np.linspace(a, b, N+1)
y = f(x)
delta_x = (b-a)/(N)

# evaluate f(x) at left, right , center points of each interval 
x_left = x[:-1]                             # left endpoints
y_left = y[:-1]
x_right = x[1:]                             # right endpoints
y_right = y[1:]
x_mid = (x[1:] + x[:-1])/2                  # mid points
y_mid = f(x_mid)


# compute approximation:
left_sum  = np.sum(y_left*delta_x)
midpt_sum = np.sum(y_mid*delta_x)
right_sum = np.sum(y_right*delta_x)
integral  = np.arctan(b) - np.arctan(a)


# function for creating textbox string
def make_txtbox_string(value_estimated, value_actual):
    textstr = '\n'.join((
    r'$\mathcal{S}=%-+.6f$' % (value_estimated, ),
    r'$\mathcal{I}=%-+.6f$' % (value_actual, ),
    r'$\mathcal{I}-\mathcal{S}=%+-.6f$' % (value_actual-value_estimated, )))
    return textstr


# plot
fig, axes= plt.subplots(1, 3, sharey=True, figsize=plt.figaspect(0.3))

# evaluation at left point
axes[0].plot(X,Y,'k')                            # function
axes[0].plot(x_left,y_left,'bo')   # numeric approximation (left value) 
axes[0].bar(x_left,y_left,width=delta_x,
       alpha=0.2,align='edge',edgecolor='b')
axes[0].axvline(x=0,color='k', linestyle=':')
axes[0].set_title('Left Sum (N = %i)'%(N))
textstr_left = make_txtbox_string(left_sum, integral)
axes[0].text(3, 0.95, textstr_left, 
             verticalalignment='top', horizontalalignment='right')
# evaluation at mid point
axes[1].plot(X,Y,'k')                            # function
axes[1].plot(x_mid,y_mid,'go')   # numeric approximation (left value) 
axes[1].bar(x_mid,y_mid,width=delta_x,
       alpha=0.2,edgecolor='g')
axes[1].axvline(x=0,color='k', linestyle=':')
axes[1].set_title('Midpoint Sum (N = %i)'%(N))
textstr_midpt = make_txtbox_string(midpt_sum, integral)
axes[1].text(3, 0.95, textstr_midpt, 
             verticalalignment='top', horizontalalignment='right')
# evaluation at right point
axes[2].plot(X,Y,'k')                            # function
axes[2].plot(x_right,y_right,'ro')   # numeric approximation (left value) 
axes[2].bar(x_right,y_right,width=-delta_x,
       alpha=0.2,align='edge',edgecolor='r')
axes[2].axvline(x=0,color='k', linestyle=':')
axes[2].set_title('Right Sum (N = %i)'%(N))
textstr_right = make_txtbox_string(right_sum, integral)
axes[2].text(3, 0.95, textstr_right, 
             verticalalignment='top', horizontalalignment='right')

plt.show()




```

##### Riemann Sums computed over Analytic Functions

---
**Exercise (1):**

Define a function that computes *left*, *right* and *midpoint* rules. The function should take 5 parameters:
-  *f*:    one-variable vectorized function
- *a*, *b*: scalar, lower and upper bound of integration interval $[a, b]$
- *N*:    scalar, number of subintervals for partitioning $[a, b]$
- *method*: string, type of sum, either 'left', 'right', 'midpoint'
   
---



```python
def riemann_sum(f, a, b, N, method='left'):
  """
  Returns riemann sum of function `f` over interval [a, b].
  
  Args:
    - f:    one-variable function
    - a, b: scalar, lower and upper bound of integration interval
    - N:    scalar, number of subintervals for partitioning [a, b]
    - method: string, type of sum, either 'left', 'right', 'midpoint'
    
  Returns:
    - integral approximation based on left, right, or midpoint method.
  """
  delta_x = np.abs(b-a)/N
  x = np.linspace(a, b, N+1)

  x_left = x[:-1]
  x_right = x[1:]
  if method == 'left':
    return np.sum(f(x_left)*delta_x)
  elif method == 'right':
    return np.sum(f(x_right)*delta_x)
  elif method == 'midpoint':
    x_midpt = (x_right + x_left)/2 
    return np.sum(f(x_midpt)*delta_x)
  else:
    print("Method '%s' is undefined."%method)
    
```

##### Convergence Behavior



---

**Exercise (2):**

Use the function defined in the previous exercise to compare the convergence behaviour for *left*, *right* and *midpoint* quadratures.

- Choose a function $f(x)$ that can be integrated analytically. Compute the definite integral $\mathcal{I}$ of that function on the interval $[a,b]$.
- Approximate the definite integral by left, right and midpoint sum using some number $N$ of partitioning intervals . Compute the approximation error.
- Repeat this procedure for different partitionings. Start with $N=1$ and double the number of intervals until $\Delta_i \approx 10^{-5}$.
- Plot the approximation error in function of $\Delta_i$, use logarithmic axes.
- Compute the convergence order using formula (7) in from the [numerical differentiation notebook](https://github.com/cohmathonc/biosci670/blob/master/IntroductionComputationalMethods/03_IntroCompMethods_NumericalDifferentiation.ipynb).
---




```python
# define function
# we use again f(x)=1/(1+x2) with antiderivative F(x)=arctan(x)
def f(x):
  return 1/(1+np.power(x,2))

# integration interval
a = 0
b = 3

# actual value of integral
integral = np.arctan(b) - np.arctan(a) 

# compute error for different partitions
# we start with N=1 and double the number of intervals, i.e.
# half the distance between partition points in each iteration
dfs = {}
for method in ['left', 'right', 'midpoint']:
  df = pd.DataFrame()
  for i in range(0,20):
      N = 1 * (2**i)
      delta_x = (b-a)/N
      integral_approx = riemann_sum(f, a, b, N, method)
      # keep record of parameters
      df.loc[i, 'method'] = method
      df.loc[i, 'N'] = N
      df.loc[i, 'delta_x'] = delta_x
      df.loc[i, 'integral_approx'] = integral_approx
  df['error'] = np.abs(integral - df['integral_approx'])
  dfs[method] = df
         

```


```python
# plot errors

fig = plt.figure(figsize=plt.figaspect(0.5))
ax = fig.add_subplot(111) 

for method in ['left', 'right', 'midpoint']:
  df = dfs[method]
  ax.loglog(df.delta_x.values, df.error.values, marker='o', 
            linestyle='-', label=method)

ax.set_xlabel('$\Delta_i$')
ax.set_ylabel("absolute error")
ax.set_title('Error in Integral Approximations using Rectangular Quadratures')
ax.legend()
plt.show()
```


```python
# we used this function previously in the numerical differentiation module

def compute_convergence_order(eps):
  """
  Computes ratio of errors steps: ratio(n) = error(n-1)/error(n)
  """
  eps_shifted = np.roll(eps, shift=1)
  eps_shifted[0] = np.nan
  eps_ratio = eps_shifted / eps
  convergence_order = np.log2(eps_ratio)
  return convergence_order
```


```python
# convergence order
for method in ['left', 'right', 'midpoint']:
  df = dfs[method]
  print("Convergence order of method '%s': "%(method), 
       compute_convergence_order(df.error.values))

```

#### Riemann Sums computed over Grids

Until now we have a assumed that an analytic form of the integrand $f(x)$ is known that can be evaluated at arbitrary points $x_i$.
This is not always the case. In fact, most often only the values $f(x_i)$ and positions $x_i$ at which  $f(x)$ has been sampled will be known, but not the function $f(x)$ itself.

---
**Exercise (3):**

1. Rewrite the numerical integration function from the previous exercise to take as input parameters 
   - an array of sampled function values $f(x_i)$
   - an array of sampling positions $x_i$
 
  Only consider *left* and *right rules*. 
 Assume that the integral is to be performed over the entire array of function values, i.e. assume that the integration bounds correspond to $x_0$ and $x_N$.

2. Why is implementing the midpoint rule not straightforward?

3. If we assume that the function $f(x)$ has been sampled on a *regular grid* with equal spacings between sampling points $\Delta_i=x_{i+1}-x_i$, then *left* and *right* sums can be computed from  $\Delta_i$ without knowing each $x_i$ in the partition.

  How could you extend your function to accept as input either an array of sampling positions $x_i$, or  the spacing between sampling points $\Delta_i$.

---




```python
def riemann_sum_grid(data, x, method='left'):
  """
  Returns riemann sum of data/function values sampled at positions x_i.
  
  Args:
    - data: an array of data points / function values
    - x:    can be a scalar, indicating the spacing between subsequent 
            observation/evaluation points, or
            an array of same length as data, indicating the observation points
    - method: string, type of sum, either 'left', or 'right'
  """    
  if isinstance(x, np.ndarray):
    delta_x = np.roll(x, -1) - x       # dx = x_i+1 - x_i
    delta_x = delta_x[:-1]             # for i = 0 ... N-1
  elif (isinstance(x, float) or isinstance(x, int)):
    delta_x = x
    
  if method == 'left':
    data_sel = data[:-1]
  elif method == 'right':
    data_sel = data[1:]
  
  riemann_sum = np.sum(delta_x * data_sel)
  return riemann_sum
```


```python
# define function to be integrated
# we use again f(x)=1/(1+x2) with antiderivative F(x)=arctan(x)
def f(x):
  return 1/(1+np.power(x,2))

# integration interval
a = 0
b = 3

# actual value of integral
integral = np.arctan(b) - np.arctan(a) 

# compute error for different partitions
# we start with N=1 and double the number of intervals, i.e.
# half the distance between partition points in each iteration
dfs = {}
for method in ['left', 'right']:
  df = pd.DataFrame()
  for i in range(0,20):
      # partition for numeric approximation
      N = 1 * (2**i)
      x = np.linspace(a, b, N+1)
      y = f(x)
      delta_x = (b-a)/(N)
      integral_approx = riemann_sum_grid(y, x, method)
      # keep record of parameters
      df.loc[i, 'method'] = method
      df.loc[i, 'N'] = N
      df.loc[i, 'delta_x'] = delta_x
      df.loc[i, 'integral_approx'] = integral_approx
  df['error'] = np.abs(integral - df['integral_approx'])
  dfs[method] = df
```


```python
# plot errors

fig = plt.figure(figsize=plt.figaspect(0.5))
ax = fig.add_subplot(111) 

for method in ['left', 'right']:
  df = dfs[method]
  ax.loglog(df.delta_x.values, df.error.values, marker='o', 
            linestyle='-', label=method)

ax.set_xlabel('$\Delta_i$')
ax.set_ylabel("absolute error")
ax.set_title('Error in Integral Approximations using Rectangular Quadratures')
ax.legend()
plt.show()
```


```python
# convergence order
for method in ['left', 'right']:
  df = dfs[method]
  print("Convergence order of method '%s': "%(method), 
       compute_convergence_order(df.error.values))
```

### Quadrature based on Interpolation

We have seen previously that the midpoint quadrature converges to the true value of the integral more quickly than the rectangular quadrature.
However, in situations when numeric integration is most useful,  the 'midpoint' value is not known nor can it be evaluated.

Since we do not know the functional form of $f(x)$ for $x\in[x_i, x_{i+1}]$, we can only approximate its values by *interpolation*.
Indeed, a whole class of quadrature rules is *based on interpolating functions*. Instead of computing the integral of $f(x)$, we may approximate the true value of the integral by integrating an interpolant of $f(x)$. 

An important class of interpolants are polynomials: 
Given the values $f(x_0), \ldots, f(x_N)$ on $N+1$ points, any function $f(x)$ can be interpolated by a polynomial $P_N(x)$ of degree $\leq N$. 

We have already used this approach:
For the left, right and midpoint *rectangular quadrature* rules, we simply approximated the function $f(x)$ in the interval $x\in[x_i, x_{i+1}]$ by a constant, i.e. a polynomial of degree 0, corresponding to the value of the function at either $x_i$, $x_{i+1}$, or the midpoint between both.


#### Trapezoidal Rule

Instead of approximating the function's value with information from only a single point $x_i$, we could use information from two neighbouring points. 
This would allow us to interpolate $f(x)$ by a polynomial of degree 1, i.e. a linear function that connects both endoints of the interval $[x_i, x_{i+1}]$. 

What would such a linear interpolation look like?

For readibility we rename the boundaries of the integration (sub)interval to $a$ and $b$, and we seek a polynomial of degree 1, $P_1(x)$, that interpolates the function f(x) for $x\in [a, b]$.

Surely, for $x=a$ and $x=b$, the interpolating function should recover the actual function values, i.e. $P_1(a)=f(a)$ and $P_1(b)=f(b)$. 
At all positions $x \in (a, b)$ for which no information about the value of $f(x)$ is available, we weight the contribution of $f(a)$ and $f(b)$ by the distance of $x$ from points $a$ and $b$ repectively:
$$f(x) \approx P_1(x) = f(a) \frac{b-x}{b-a} + f(b)\frac{x-a}{b-a} \qquad x\in[a,b] \tag{4}$$

Having found an interpolation of $f(x)$, we can now approximate the integral:

\begin{align}
\int_a^b f(x)\, dx \approx \int_a^b P_1(x)\, dx &= f(a) \int_a^b \frac{b-x}{b-a} \, dx+ f(b) \int_a^b \frac{x-a}{b-a}\, dx \\
                                                                                &=  f(a) \frac{b-a}{2} + f(b)\frac{b-a}{2} \\
                                                                                &= \frac{b-a}{2}\left(f(a) + f(b)\right) \tag{5}
\end{align}

Geometrically, this means that the 'area under the curve' in each subinterval $[a,b]=[x_i, x_{i+1}]$ is approximated, not by a rectangle as before, but by a trapezoid.
This quadrature is therefore called the **trapezoidal rule** or **trapezoidal quadrature**.

Applied to the partition $\{x_0, x_1, \ldots x_{N}\}$ from before:

\begin{align}
\text{trapezoid sum} \qquad \mathcal{S}_{T, N}  &=\sum_{i=0}^{N-1} \frac{x_{i+1} - x_i}{2}\, \left[ f(x_{i+1}) + f(x_{i}) \right] \tag{6a}\\
                                                                                    &= \frac{x_{1} - x_0}{2}\, \left[ f(x_{1}) + f(x_{0}) \right] 
                                                                                      + \frac{x_{2} - x_1}{2}\, \left[ f(x_{2}) + f(x_{1}) \right] 
                                                                                      + \dots
                                                                                      + \frac{x_{N-1} - x_{N-2}}{2}\, \left[ f(x_{N-1}) + f(x_{N-2}) \right]
                                                                                      + \frac{x_{N} - x_{N-1}}{2}\, \left[ f(x_{N}) + f(x_{N-1}) \right]\\
                                                                                   &= f(x_{0})\frac{\Delta x}{2} +  2 f(x_{1})\frac{\Delta x}{2} + 
                                                                                        2 f(x_{2})\frac{\Delta x}{2} +\dots + 2 f(x_{N-2})\frac{\Delta x}{2} + 
                                                                                        2 f(x_{N-1})\frac{\Delta x}{2} + f(x_{N})\frac{\Delta x}{2}\\
                                                                                   &=  f(x_{0})\frac{\Delta x}{2} 
                                                                                       +  \sum_{i=1}^{N-1} f(x_i)\, \Delta x 
                                                                                       + f(x_{N})\frac{\Delta x}{2} \tag{6b}
\end{align}





---
**Exercise (4):**

1. Define a function that uses the *trapezoidal quadrature* to compute an integral approximation.
  The function should take 2 parameters:

   - *data*: an array of sampled function values $f(x_i)$
   - *x*: an array of sampling positions $x_i$
 
  Assume that the integral is to be performed over the entire array of function values, i.e. assume that the integration bounds correspond to $x_0$ and $x_N$.

2.  Extend your function to accept as input either an array of sampling positions $x_i$, or  the spacing between sampling points $\Delta_i$, i.e. to compute (6a) and (6b).

---





```python
def trapezoidal_sum_grid(data, x):
  """
  Returns riemann sum of data / function using trapezoidal quadrature.
  
  Args:
    - data: an array of data points / function values
    - x:    can be a scalar, indicating the spacing between subsequent 
            observation/evaluation points, or
            an array of same length as data, indicating the observation points
    - method: string, type of sum, either 'left', or 'right'
  """  
  if isinstance(x, np.ndarray):
    data_i_plus_iplus1 = data + np.roll(data,-1)  # f(x_i) + f(x_i+1)
    data_i_plus_iplus1 = data_i_plus_iplus1[:-1]  # for i = 0 ... N-1
    delta_x = np.roll(x, -1) - x                  # dx = x_i+1 - x_i
    delta_x = delta_x[:-1]                        # for i = 0 ... N-1
    return np.sum(0.5 * data_i_plus_iplus1 * delta_x)
  elif (isinstance(x, float) or isinstance(x, int)):
    delta_x = x
    weights = np.ones_like( data )
    weights[0] = 0.5
    weights[-1]= 0.5
    return np.sum(delta_x * weights * data)

```


```python
# define function
# we use again f(x)=1/(1+x2) with antiderivative F(x)=arctan(x)
def f(x):
  return 1/(1+np.power(x,2))

# integration interval
a = 0
b = 3

# actual value of integral
integral = np.arctan(b) - np.arctan(a) 

# compute error for different partitions
# we start with N=1 and double the number of intervals, i.e.
# half the distance between partition points in each iteration
df = pd.DataFrame()
for i in range(0,20):
    # partition for numeric approximation
    N = 1 * (2**i)
    x = np.linspace(a, b, N+1)
    y = f(x)
    delta_x = (b-a)/(N)
    integral_approx = trapezoidal_sum_grid(y, x)
    #integral_approx = trapezoidal_sum_grid(y, delta_x)
    
    # keep record of parameters
    df.loc[i, 'N'] = N
    df.loc[i, 'delta_x'] = delta_x
    df.loc[i, 'integral_approx'] = integral_approx
df['error'] = np.abs(integral - df['integral_approx'])

```


```python
# plot errors

fig = plt.figure(figsize=plt.figaspect(0.5))
ax = fig.add_subplot(111) 

ax.loglog(df.delta_x.values, df.error.values, marker='o', 
          linestyle='-', label='trapezoidal quadrature')

ax.set_xlabel('$\Delta_i$')
ax.set_ylabel("absolute error")
ax.set_title('Error in Integral Approximations using Trapezoidal Quadrature')
ax.legend()
plt.show()
```


```python
# convergence order
print("Convergence order of trapezoidal quadrature: ", 
     compute_convergence_order(df.error.values))
```

#### Higher Order Quadratures

The procedure we have used to derive the trapezoidal rule above also allows to derive higher order quadratures by first finding an interpolant that appoximates $f(x)$ on the interval $[a, b]$ and then integrating that interpolant in place of $f(x)$.
Standard methods exist for finding interpolants of a function, given its values at specific locations. The approach demonstrated for the derivation of the trapezoidal used a first degree [Lagrange Polynomial](https://en.wikipedia.org/wiki/Lagrange_polynomial).

## Useful Python / Numpy Functions


```python
# np.trapz implements the trapezoidal rule

import numpy as np
# we use again f(x)=1/(1+x2) with antiderivative F(x)=arctan(x)
def f(x):
  return 1/(1+np.power(x,2))

# integration interval
a = 0
b = 3

# partition
N = 100
delta_x = (b-a)/(N)
x = np.linspace(a, b, N+1)

# values
y = f(x)

# integration using x
integral_trapez_our = trapezoidal_sum_grid(y, x)
integral_trapez_np  = np.trapz(y, x=x)

print("our implementation (given array x): ", integral_trapez_our)
print("np.trapz implementation (given array x): ", integral_trapez_np)

# integration using delta_x
integral_trapez_our = trapezoidal_sum_grid(y, delta_x)
integral_trapez_np  = np.trapz(y, x=None, dx=delta_x)

print("our implementation (given delta x): ", integral_trapez_our)
print("np.trapz implementation (given delta x): ", integral_trapez_np)
```


```python
# scipy.integrate: trapz, simps
import scipy.integrate as integrate


import numpy as np
# we use again f(x)=1/(1+x2) with antiderivative F(x)=arctan(x)
def f(x):
  return 1/(1+np.power(x,2))

# integration interval
a = 0
b = 3

# partition
N = 100
delta_x = (b-a)/(N)
x = np.linspace(a, b, N+1)

# values
y = f(x)

# integration using x
integral_trapez_our = trapezoidal_sum_grid(y, x)
integral_trapz = integrate.trapz(y, x)
#integral_cumtrapz = integrate.cumtrapz(y, x)
integral_simps = integrate.simps(y, x)

print("our implementation: ", integral_trapez_our)
print("trapz   : ", integral_trapz)     # uses linear interpolant as discussed
#print("cumtrapz: ", integral_cumtrapz)
print("simps   : ", integral_simps)     # uses quadratic interpolant


# above functions operate on discrete observations
# scipy.integrate also provides functions for numeric integration provided a 
# function object
```

## Exercises

- In [this](https://github.com/cohmathonc/biosci670/blob/master/IntroductionComputationalMethods/exercises/05_NumericalIntegration.ipynb) exercise you approximate an definite integral using the trapezoidal quadrature. (**Optional**) 
- In [this](https://github.com/cohmathonc/biosci670/blob/master/IntroductionComputationalMethods/exercises/06_NumericalIntegration_ComputeErrorFunction.ipynb) exercise you compute the *error function* using numerical integration. (**Optional**) 

###### About 
This notebook is part of the *biosci670* course on *Mathematical Modeling and Methods for Biomedical Science*.
See https://github.com/cohmathonc/biosci670 for more information and material.

7628fac794ced560c928c63f034cd46fca5de7c7
ipynb
Jupyter Notebook
IntroductionComputationalMethods/04_IntroCompMethods_NumericalIntegration.ipynb
sbranciamore/biosci670coh
92e91aa1e0266b54b0831537b857796eaf0b949c
CC-BY-4.0
2019-02-19T21:20:42.000Z
2021-01-12T20:25:42.000Z
IntroductionComputationalMethods/04_IntroCompMethods_NumericalIntegration.ipynb
danielabler/biosci670
a9fc3e3062cc7ee4f2da3c26a4c5fc22c5cb8cb4
CC-BY-4.0
IntroductionComputationalMethods/04_IntroCompMethods_NumericalIntegration.ipynb
danielabler/biosci670
a9fc3e3062cc7ee4f2da3c26a4c5fc22c5cb8cb4
CC-BY-4.0
2019-05-15T17:33:41.000Z
2021-02-22T18:11:37.000Z
Qwen/Qwen-72B
1. YES
2. YES