```python
%pylab inline
import numpy as np
import pandas as pd
import sympy as sp
```

    Populating the interactive namespace from numpy and matplotlib



```python
font = {'size'   : 14}
matplotlib.rc('font', **font)
sp.init_printing()
```

# Approximate integration methods

Recall that the definite integral is defined as a limit of Riemann sums, so any Riemann sum could be used as an approximation to the integral: If we divide $[a,b]$ into $n$ subintervals of equal length $\Delta x = (b-a)/n$, then we have

$$ \int_a^b f(x)\ dx \approx \sum\limits_{i=1}^n f(x_i^*) \Delta x$$

where $x_i^*$ is any point in the $i$th subinterval $[x_{i-1}, x_i]$. 

## Left endpoint approximation

If $x_i^*$ is chosen to be the left-endpoint of the interval, then 

$$ \tag{1} \int_a^b f(x)\ dx \approx L_n = \sum\limits_{i=1}^n f(x_{i-1}^*) \Delta x$$

This is the implementation of equation $(1)$.


```python
def L(f, a, b, n):
    dx = (b - a) / n
    X = np.linspace(a, b - dx, n)
    return sum([f(x)*dx for x in X])
```

## Right endpoint approximation

If we choose $x_i^*$ to be the right endpoint, then $x_i^*=x_i$ and we have

$$ \tag{2} \int_a^b f(x)\ dx \approx R_n = \sum\limits_{i=1}^n f(x_{i}^*) \Delta x$$

The approximations $L_n$ and $R_n$ are called the **left endpoint approximation** and **right endpoint approximation** respectively. Next is the implementation of equation $(2)$.


```python
def R(f, a, b, n):
    dx = (b - a) / n
    X = np.linspace(a, b - dx, n) + dx
    return sum([f(x)*dx for x in X])
```

## Midpoint rule

Let's look at the midpoint approximation $M_n$, which appears to be better than either $L_n$ or $R_n$.

$$ \tag{3} \int_a^b f(x)\ dx \approx M_n = \Delta x[f(\bar{x}_1), f(\bar{x}_2), \ldots, f(\bar{x}_n)] $$ where $$\Delta x = \frac{b-a}{n}$$ and $$\bar{x}_i = \frac{1}{2}(x_{i-1} + x_i) $$ is the midpoint of $[x_{i-1}, x_i]$. Equation $(3)$ also has an implementation.


```python
def M(f, a, b, n):
    dx = (b - a) / n
    X = np.linspace(a, b - dx, n) + dx / 2
    return sum([f(x)*dx for x in X])
```

## Trapezoidal rule

Another approximation, called the Trapezoidal Rule, results from averaging the appromations in equations 1 and 2.

$$ \int_a^b f(x)\ dx \approx \dfrac{1}{2}\left[ \sum\limits_{i=1}^n f(x_{i-1})\Delta x + \sum\limits_{i=1}^n f(x_{i})\Delta x \right] = \dfrac{\Delta x}{2} \sum\limits_{i=1}^n \left( f(x_{i-1})+f(x_{i})\right) $$


```python
def T(f, a, b, n):
    dx = (b - a) / n
    X = np.linspace(a, b - dx, n)
    return sum( [f(x) + f(x+dx) for x in X] ) * dx / 2
```

It can also be defined as:

$$ \tag{4} \int_a^b f(x)\ dx \approx T_n = \dfrac{\Delta x}{2} \left[ f(x_0) + 2f(x_1) + 2f(x_2) + \ldots + 2f(x_{n-1}) + f(x_n)\right]$$

where $\Delta x = (b-a)/n$ and $x_i=a+i\Delta x$.

# Approximations

Suppose we have the following integral:

$$\tag{Example 1}\int^2_1 \dfrac{1}{x}\ dx = \begin{bmatrix} \ln\ |\ x\ | \end{bmatrix}^2_1 = \ln 2 - \ln 1 = \ln 2 \approx 0.693147.$$

It is an easy integral, because we want to check the error between the approximations methods and the exact value. Estimating the integral with each of the methods gives:


```python
f=lambda x: 1/x
a=1; b=2
N=[]; Ln=[]; Rn=[]; Mn=[]; Tn=[];

for n in [5, 10, 15, 20, 25]:
    N.append(n)
    Ln.append(L(f, a, b, n))
    Rn.append(R(f, a, b, n))
    Mn.append(M(f, a, b, n))
    Tn.append(T(f, a, b, n))
```

Now we add everything to a datatable and calculate the error between the exact value of the integral and each method.


```python
df=pd.DataFrame()
df['n'] = N
df['ln 2'] = [math.log(2)] * len(N)
df['Ln'] = Ln
df['Rn'] = Rn
df['Mn'] = Mn
df['Tn'] = Tn
df['L_error'] = math.log(2) - df['Ln'] 
df['R_error'] = math.log(2) - df['Rn']
df['M_error'] = math.log(2) - df['Mn']
df['T_error'] = math.log(2) - df['Tn']
df
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>n</th>
      <th>ln 2</th>
      <th>Ln</th>
      <th>Rn</th>
      <th>Mn</th>
      <th>Tn</th>
      <th>L_error</th>
      <th>R_error</th>
      <th>M_error</th>
      <th>T_error</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>5</td>
      <td>0.693147</td>
      <td>0.745635</td>
      <td>0.645635</td>
      <td>0.691908</td>
      <td>0.695635</td>
      <td>-0.052488</td>
      <td>0.047512</td>
      <td>0.001239</td>
      <td>-0.002488</td>
    </tr>
    <tr>
      <th>1</th>
      <td>10</td>
      <td>0.693147</td>
      <td>0.718771</td>
      <td>0.668771</td>
      <td>0.692835</td>
      <td>0.693771</td>
      <td>-0.025624</td>
      <td>0.024376</td>
      <td>0.000312</td>
      <td>-0.000624</td>
    </tr>
    <tr>
      <th>2</th>
      <td>15</td>
      <td>0.693147</td>
      <td>0.710091</td>
      <td>0.676758</td>
      <td>0.693008</td>
      <td>0.693425</td>
      <td>-0.016944</td>
      <td>0.016389</td>
      <td>0.000139</td>
      <td>-0.000278</td>
    </tr>
    <tr>
      <th>3</th>
      <td>20</td>
      <td>0.693147</td>
      <td>0.705803</td>
      <td>0.680803</td>
      <td>0.693069</td>
      <td>0.693303</td>
      <td>-0.012656</td>
      <td>0.012344</td>
      <td>0.000078</td>
      <td>-0.000156</td>
    </tr>
    <tr>
      <th>4</th>
      <td>25</td>
      <td>0.693147</td>
      <td>0.703247</td>
      <td>0.683247</td>
      <td>0.693097</td>
      <td>0.693247</td>
      <td>-0.010100</td>
      <td>0.009900</td>
      <td>0.000050</td>
      <td>-0.000100</td>
    </tr>
  </tbody>
</table>
</div>




```python
f, ((ax1, ax2, ax3, ax4)) = plt.subplots(1,4, sharex='col', sharey='row'
                                         , figsize=(14,6))
cols = ['L_error', 'R_error', 'M_error', 'T_error']
titles = ['Left endpoint', 'Right endpoint', 'Midpoint', 'Trapezoidal']
axes = [ax1, ax2, ax3, ax4]

for i in range(len(cols)):
    ax = axes[i]
    ax.bar(df['n'], df[cols[i]], width=3, alpha=0.33, fc='b')
    ax.set_title(titles[i])
    ax.set_ylabel('Error')
    ax.set_xlabel('$n$ approximations')
    ax.axhline(0, c='black', lw=1, ls='dashed')
    ax.grid(ls='dashed', alpha=0.5)
```

We can make several observations from these tables:

1. In all methods we get more accurate approximations when we increase the value of $n$. However, very large numbers of $n$ result in so many arithmetic operations that we have to beware of accumulated round-off error. 
2. The errors in the left and right endpoint approximations are opposite in sign and appear to decrease by a factor of about 2 when we double the value of $n$.
3. The Trapezoidal and Midpoint Rules are much more accurate that the endpoint approximations.
4. The errors in the Trapezoidal and Midpoint Rules are opposite in sign and appear to decrease by a factor of about 4 when we double the value of $n$.
5. The size of error in the Midpoint Rule is about half the size of the error in the Trapezoidal Rule.

# Error bounds

Suppose $|\ f''(x)\ | \leq K$ for $a \leq x \leq b$. If $E_t$ and $E_m$ are the errors in the Trapezoidal and Midpoint Rules, then $$\tag{5} |\ E_t\ | \leq \dfrac{K(b-a)^3}{12n^2} \quad \text{and} \quad |\ E_m\ | \leq \dfrac{K(b-a)^3}{24n^2}.$$ Let's apply this error estimate to the Trapezoidal Rule approximation in Example 1 equation. First we find $f''(x)$:


```python
x = sp.symbols('x')
f = 1/x
sp.diff(sp.diff(f))
```

Because $1 \leq x \leq 2$, we have $1/x \leq 1$, so

$$ |\ f''(x)\ | = \left|\dfrac{2}{x^3}\right| \leq \dfrac{2}{1^3}=2 $$

Therefore, taking $K=2,a=1,b=2$, and $n=5$ in the error estimate equation 5, we see that

$$ |\ E_r |\ \leq \dfrac{2(2-1)^3}{12(5)^2}=\dfrac{1}{150}\approx0.006667 $$

**Guaranteed error** How large should we take $n$ in order to guarantee that the Trapezoidal and Midpoint Rule approximations for $\int_1^2(1/x)\ dx$ are accurate to within $0.0001$?

We saw in the preceding calculation that $|\ f(x)\ | \leq 2$ for $1\leq x\leq 2$, so we take $K=2, a=1$, and $b=2$ in equation 5. Accuracy to within $0.0001$ means that the size of the error should be less than $0.0001$. Therefore we choose $n$ so that $$ \dfrac{2(1)^3}{12n^2}\leq0.0001.$$

Solving the inequality for $n$, we get $$ n \gt \dfrac{1}{\sqrt{0.0006}}\approx 40.8 $$ Thus $n=41$ will ensure the desired accuracy.

# Simpson's rule

...


```python
from sympy.integrals.manualintegrate import integral_steps
```


```python
x = sp.symbols('x')
integral_steps(x*sp.sin(3*x), x)
```




    PartsRule(u=x, dv=sin(3*x), v_step=URule(u_var=_u, u_func=3*x, constant=1/3, substep=ConstantTimesRule(constant=1/3, other=sin(_u), substep=TrigRule(func='sin', arg=_u, context=sin(_u), symbol=_u), context=sin(_u), symbol=_u), context=sin(3*x), symbol=x), second_step=ConstantTimesRule(constant=-1/3, other=cos(3*x), substep=URule(u_var=_u, u_func=3*x, constant=1/3, substep=ConstantTimesRule(constant=1/3, other=cos(_u), substep=TrigRule(func='cos', arg=_u, context=cos(_u), symbol=_u), context=cos(_u), symbol=_u), context=cos(3*x), symbol=x), context=-cos(3*x)/3, symbol=x), context=x*sin(3*x), symbol=x)




```python

```

f0f781cff1fa356c39a06b304720557e8a48469b
ipynb
Jupyter Notebook
Notebooks/Approximate integration.ipynb
darkeclipz/jupyter-notebooks
5de784244ad9db12cfacbbec3053b11f10456d7e
Unlicense
2018-08-28T12:16:12.000Z
2018-08-28T12:16:12.000Z
Notebooks/Approximate integration.ipynb
darkeclipz/jupyter-notebooks
5de784244ad9db12cfacbbec3053b11f10456d7e
Unlicense
Notebooks/Approximate integration.ipynb
darkeclipz/jupyter-notebooks
5de784244ad9db12cfacbbec3053b11f10456d7e
Unlicense
Qwen/Qwen-72B
1. YES
2. YES