Anthony Jimenez
25 August 2021
http://www.fa-jimenez.com/
import pandas as pd
import numpy as np
from scipy.integrate import quad
import matplotlib.pyplot as plt
%matplotlib inline
This notebook will cover how one can manually program a numerical integration function that leverages rectangular panels to integrate any given function.
The main libraries that will be used are numpy (for arrays) and scipy for pre-built integration functions (for comparison purposes).
Given an $x$ domain and $n$ panels desired from the input, the given function will be integrated $F(x) = \int_{a}^{b} f(x)\, dx$.
In this notebook we wil look at the domain, $x: [0, 9]$ for the following function:
$$f(x) = \text{sin}(x)$$$$F(x) = \int_{a}^{b} f(x)\, dx = -\text{cos}(x)|_{a}^{b}$$When evaluating this integral over the aforementioned domain, we get the following result.
$$F(x) = \int_{0}^{9} \text{sin}(x)\, dx = -\text{cos}(x)|_{0}^{9} = 1.91113026188468$$# Defining the domain and function to map
x = np.linspace(0, 9, 1000)
f = lambda x: np.sin(x)
y = list(map(f, x))
# Plot the sin() function
fig, ax = plt.subplots(figsize=[12,4])
ax.plot(x, y)
ax.minorticks_on()
ax.set_aspect('equal')
ax.set_xlabel('x')
ax.set_ylabel('y')
plt.title(r'Domain of $f(x) = $sin$(x)$ to be investigated')
Text(0.5, 1.0, 'Domain of $f(x) = $sin$(x)$ to be investigated')
Provided the function, number of panels the rectangular method is computed as follows:
$$w_{\text{panel}} = \frac{b - a}{n}$$$$h = f(x)$$$$A_{\text{rectangle}} = h * w_{\text{panel}}$$The main difference between the two rectangular methods will be which $x$ values will be used for evaluating the height of the rectangle, $h$. When implementing the left rectangular method, the height is always evaluated at the left-hand point of any given rectangle panel. As a result, the last $x$ value will not be used when computing the height of the rectangle.
In contrast the right rectangular method will skip the first $x$ value of the domain and the second value will be used to compute the height of the first rectangle.
def rectangular_integration(f, n_panels, x_interval, left_or_right):
# Compute the panel start/end points and rectangular width
xpts_all = np.linspace(start=x_interval[0], stop=x_interval[1], num=n_panels+1)
width = xpts_all[1] - xpts_all[0]
# Determine if using left rectangle point or right
if left_or_right == 'left':
xpts = xpts_all[:-1]
elif left_or_right == 'right':
xpts = xpts_all[1:]
# Compute the integral function
ypts = np.array(list(map(f, xpts)))
fpts = width * ypts
Fpts = np.cumsum(fpts)
integral_output = Fpts[-1]
return integral_output, xpts, ypts, fpts, xpts_all, Fpts
Naturally, with lower density rectangle panels, each rectangle has a larger width. This means that the height computed for the rectangle will be applied to a larger domain section and this can lead to bad estimates if the function changes quickly.
In effect, a lower density count of rectangles will lead to the function being averaged over a bigger section. Below the true function is shown with subplots also showing the interpretation of the rectangles using both left and right schemes.
Lastly, a print out of the integral result for left- and right-hand side are displayed for comparison.
# Set the domain to integrate over and number of rectangles
x_domain = [0, 9]
n_panels = 11
F_left = rectangular_integration(f, n_panels, x_domain, 'left')
F_right = rectangular_integration(f, n_panels, x_domain, 'right')
print(f'Left rectangular integration: {F_left[0]}')
print(f'Right rectangular integration: {F_right[0]}')
fig, axs = plt.subplots(3, figsize=[20,15])
axs[0].plot(x, y, label='Function')
axs[1].plot(x, y, label='Function')
axs[1].bar(F_left[1], F_left[2], linewidth=0.01, color='orange', alpha=0.5, label='Left Rectangle')
axs[2].plot(x, y, label='Function')
axs[2].bar(F_right[1], F_right[2], linewidth=0.01, color='lightcoral', alpha=0.5, label='Right Rectangle')
axs[0].set_xlim(x_domain[0], x_domain[1])
axs[1].set_xlim(x_domain[0], x_domain[1])
axs[2].set_xlim(x_domain[0], x_domain[1])
axs[0].legend()
axs[1].legend()
axs[2].legend()
Left rectangular integration: 1.6347149362849438 Right rectangular integration: 1.971902787846381
<matplotlib.legend.Legend at 0x26812ef1dc0>
With increased rectangles, we see small width panels and the height is thus averaged over a smaller section. When comparing this with the true function, we see that the shape of the sin() function is a lot more alike and we can see that the integral approximations are getting closer to the analytical answer.
# Set the domain to integrate over and number of rectangles
x_domain = [0, 9]
n_panels = 30
F_left = rectangular_integration(f, n_panels, x_domain, 'left')
F_right = rectangular_integration(f, n_panels, x_domain, 'right')
print(f'Left rectangular integration: {F_left[0]}')
print(f'Right rectangular integration: {F_right[0]}')
fig, axs = plt.subplots(3, figsize=[20,15])
axs[0].plot(x, y, label='Function')
axs[1].plot(x, y, label='Function')
axs[1].bar(F_left[1], F_left[2], linewidth=0.01, color='orange', alpha=0.5, label='Left Rectangle')
axs[2].plot(x, y, label='Function')
axs[2].bar(F_right[1], F_right[2], linewidth=0.01, color='lightcoral', alpha=0.5, label='Right Rectangle')
axs[0].set_xlim(x_domain[0], x_domain[1])
axs[1].set_xlim(x_domain[0], x_domain[1])
axs[2].set_xlim(x_domain[0], x_domain[1])
axs[0].legend()
axs[1].legend()
axs[2].legend()
Left rectangular integration: 1.8349574657430443 Right rectangular integration: 1.9585930113155712
<matplotlib.legend.Legend at 0x26813067550>
# Set the domain to integrate over and number of rectangles
x_domain = [0, 9]
n_panels = 1000
F_left = rectangular_integration(f, n_panels, x_domain, 'left')
F_right = rectangular_integration(f, n_panels, x_domain, 'right')
print(f'Left rectangular integration: {F_left[0]}')
print(f'Right rectangular integration: {F_right[0]}')
fig, axs = plt.subplots(3, figsize=[20,15])
axs[0].plot(x, y, label='Function')
axs[1].plot(x, y, label='Function')
axs[1].bar(F_left[1], F_left[2], width=0.03, linewidth=0.01, color='orange', alpha=0.1, label='Left Rectangle')
axs[2].plot(x, y, label='Function')
axs[2].bar(F_right[1], F_right[2], width=0.03, linewidth=0.01, color='lightcoral', alpha=0.1, label='Right Rectangle')
axs[0].set_xlim(x_domain[0], x_domain[1])
axs[1].set_xlim(x_domain[0], x_domain[1])
axs[2].set_xlim(x_domain[0], x_domain[1])
axs[0].legend()
axs[1].legend()
axs[2].legend()
Left rectangular integration: 1.909262828554406 Right rectangular integration: 1.9129718949215817
<matplotlib.legend.Legend at 0x268149b0820>
This shows the plot of the integral results for all values from [0,9]. We can see that whenever the derivative of the function approaches 0 the integration is doing the best because it is most like the shape of the rectangle. Whenever the function is growing or declining that is when we see the integration missing in a worse manner.
# Recall domain of the problem we are interested in
x_domain = [0, 9]
x_domain_inter = np.linspace(0, 9, 30)
# Compute scipy numerical integration and actual analytical integration
F_scipy = list(map(lambda xx: quad(f, x_domain[0], xx)[0], x_domain_inter))
F_analytical = - np.cos(x_domain_inter) - (-np.cos(x_domain[0]))
# Defined the density levels to be analyzed
n_panels_low = 11
n_panels_mid = 30
n_panels_high = 1000
# Compute the left/right rectangular intergration for varying densities
F_left_low = rectangular_integration(f, n_panels_low, x_domain, 'left')
F_left_mid = rectangular_integration(f, n_panels_mid, x_domain, 'left')
F_left_high = rectangular_integration(f, n_panels_high, x_domain, 'left')
F_right_low = rectangular_integration(f, n_panels_low, x_domain, 'right')
F_right_mid = rectangular_integration(f, n_panels_mid, x_domain, 'right')
F_right_high = rectangular_integration(f, n_panels_high, x_domain, 'right')
# Plot the integral results (i.e., area underneath the sin(x) curve)
fig, axs = plt.subplots(6, figsize=[20,20])
axs[0].plot(x_domain_inter, F_analytical, label='Analytical Solution')
axs[0].scatter(x_domain_inter, F_scipy, s=120, facecolors='None', edgecolors='g', linewidths=1.7, label='Scipy Solution')
axs[0].bar(F_left_low[1], F_left_low[5], linewidth=0.01, color='orange', alpha=0.5, label='Left Rect. Low Density')
axs[1].scatter(x_domain_inter, F_scipy, s=120, facecolors='None', edgecolors='g', linewidths=1.7, label='Scipy Solution')
axs[1].plot(x_domain_inter, F_analytical, label='Analytical Solution')
axs[1].bar(F_left_mid[1], F_left_mid[5], linewidth=0.01, color='pink', alpha=0.5, label='Left Rect. Mid Density')
axs[2].plot(x_domain_inter, F_analytical, label='Analytical Solution')
axs[2].scatter(x_domain_inter, F_scipy, s=120, facecolors='None', edgecolors='g', linewidths=1.7, label='Scipy Solution')
axs[2].bar(F_left_high[1], F_left_high[5], width=0.03, linewidth=0.01, color='chocolate', alpha=0.1, label='Left Rect. High Density')
axs[3].plot(x_domain_inter, F_analytical, label='Analytical Solution')
axs[3].scatter(x_domain_inter, F_scipy, s=120, facecolors='None', edgecolors='g', linewidths=1.7, label='Scipy Solution')
axs[3].bar(F_right_low[1], F_right_low[5], linewidth=0.01, color='orange', alpha=0.5, label='Right Rect. Low Density')
axs[4].scatter(x_domain_inter, F_scipy, s=120, facecolors='None', edgecolors='g', linewidths=1.7, label='Scipy Solution')
axs[4].plot(x_domain_inter, F_analytical, label='Analytical Solution')
axs[4].bar(F_right_mid[1], F_right_mid[5], linewidth=0.01, color='pink', alpha=0.5, label='Right Rect. Mid Density')
axs[5].plot(x_domain_inter, F_analytical, label='Analytical Solution')
axs[5].scatter(x_domain_inter, F_scipy, s=120, facecolors='None', edgecolors='g', linewidths=1.7, label='Scipy Solution')
axs[5].bar(F_right_high[1], F_right_high[5], width=0.03, linewidth=0.01, color='chocolate', alpha=0.1, label='Right Rect. High Density')
axs[0].set_xlim(x_domain[0], x_domain[1])
axs[1].set_xlim(x_domain[0], x_domain[1])
axs[2].set_xlim(x_domain[0], x_domain[1])
axs[3].set_xlim(x_domain[0], x_domain[1])
axs[4].set_xlim(x_domain[0], x_domain[1])
axs[5].set_xlim(x_domain[0], x_domain[1])
axs[0].legend()
axs[1].legend()
axs[2].legend()
axs[3].legend()
axs[4].legend()
axs[5].legend()
<matplotlib.legend.Legend at 0x2681706e820>
Here we review the answers for the integration under question. We compute the percent error for each method based to the analytical solution mentioned earlier.
The quadrature method provided by the Scipy function performs the best with a percent error less than 1E-13. Due to the nature of this function (i.e., its shape) the right hand rectangular method performs better than the left hand rectangular method.
In future work, we can look at using other methods like trapezoidal integration instead of rectangles, and finally the quadrature method.
$$\text{Percent Error} = \frac{|\text{Measured Value}-\text{True Value}|}{\text{True Value}}*100\text{%}$$# Compute percent error
err_scipy = abs(F_scipy[-1] - F_analytical[-1]) / F_analytical[-1] * 100
err_left_low = abs(F_left_low[0] - F_analytical[-1]) / F_analytical[-1] * 100
err_left_mid = abs(F_left_mid[0] - F_analytical[-1]) / F_analytical[-1] * 100
err_left_high = abs(F_left_high[0] - F_analytical[-1]) / F_analytical[-1] * 100
err_right_low = abs(F_right_low[0] - F_analytical[-1]) / F_analytical[-1] * 100
err_right_mid = abs(F_right_mid[0] - F_analytical[-1]) / F_analytical[-1] * 100
err_right_high = abs(F_right_high[0] - F_analytical[-1]) / F_analytical[-1] * 100
fig, ax = plt.subplots(figsize=[20,8])
ax.bar('Low Density Left Rect.', err_left_low)
ax.bar('Mid Density Left Rect.', err_left_mid)
ax.bar('High Density Left Rect.', err_left_high)
ax.bar('Low Density Right Rect.', err_right_low)
ax.bar('Mid Density Right Rect.', err_right_mid)
ax.bar('High Density Right Rect.', err_right_high)
ax.bar('Scipy Quadrature', err_scipy)
ax.set_xlabel('\n\nIntegration Method')
ax.set_ylabel('Percent Error (%)')
plt.title('Comparison of the analytical solution to numerical integration methods')
y_plot = [err_left_low, err_left_mid, err_left_high, err_right_low, err_right_mid, err_right_high, err_scipy]
for i, v in enumerate(y_plot):
ax.text(i-0.15, v+.2, f'{v:.3f}%', color='black', fontweight='bold')
