{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"nbpresent": {
"id": "db86c5af-b769-4d51-afe0-c5b27884bf3b"
}
},
"source": [
"# Numerical Integration\n",
"## CH EN 2450 - Numerical Methods\n",
"**Prof. Tony Saad (www.tsaad.net)
Department of Chemical Engineering
University of Utah**\n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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}$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Newton-Cotes Integration\n",
"The general idea behind Newton-Cotes Integration is to approximate the function using a polynomial and then integrating the polynomial."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Left, Right, and Midpoint Rules\n",
"The left, right, and midpoint rules are the simplest of all numerical integration methods. They are based on approximating the function as a constant\n",
"\\begin{equation}\n",
"\\int_a^b f(x) \\text{d}x \\approx \\int_a^b C \\text{d}x = C\\times(b-a)\n",
"\\end{equation}\n",
"The constant can be chosen in several ways, the most common are\n",
"1. Left point rule: choose $C = f(a)$\n",
"2. Right point rule: choose $C = f(b)$\n",
"3. Midpoint rule: choose $C = f(\\frac{a+b}{2})$\n",
"\n",
"The three options are shown on the plots below."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"%matplotlib inline\n",
"%config InlineBackend.figure_format = 'svg'\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"def myf(x):\n",
" return 1.0/np.sqrt(np.pi) * np.exp(-x**2)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0.0, 0.5923990627251441)"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"xalot = np.linspace(0,2,200)\n",
"\n",
"npts = 5\n",
"\n",
"x = np.linspace(0,2, npts)\n",
"fig, (ax1, ax2, ax3) = plt.subplots(1, 3, sharey=True)\n",
"fig.set_size_inches((12,3))\n",
"a = 0.5\n",
"b = 1.0\n",
"\n",
"ax = ax1\n",
"ax.plot(xalot, myf(xalot))\n",
"ax.set_title('Leftpoint Rule')\n",
"# ax1.fill_between([a,b],[myf(a),myf(a)], edgecolor='k', facecolor='lightgray',linewidth=1)\n",
"# plot left Riemann Sum\n",
"for i in range(0, len(x)-1):\n",
" a = x[i]\n",
" b = x[i+1]\n",
" ax.fill_between([a,b],[myf(a),myf(a)], edgecolor='k', facecolor='lightgray',linewidth=1)\n",
"\n",
"ax = ax2\n",
"ax.plot(xalot, myf(xalot))\n",
"# ax3.fill_between([a,b],[myf( (a+b)/2) ,myf( (a+b)/2)], edgecolor='k', facecolor='lightgray',linewidth=1)\n",
"# plot midpoint rule\n",
"ax.set_title('Midpoint Rule')\n",
"for i in range(0, len(x)-1):\n",
" a = x[i]\n",
" b = x[i+1]\n",
" h = myf((a+b)/2)\n",
" ax.fill_between([a,b],[h,h],edgecolor='k', facecolor='lightgray',linewidth=1)\n",
"\n",
"ax = ax3\n",
"ax.plot(xalot, myf(xalot))\n",
"# ax2.fill_between([a,b],[myf(b),myf(b)], edgecolor='k', facecolor='lightgray',linewidth=1)\n",
"# plot right Riemann Sum\n",
"ax.set_title('Rightpoint Rule')\n",
"for i in range(0, len(x)-1):\n",
" a = x[i]\n",
" b = x[i+1]\n",
" ax.fill_between([a,b],[myf(b),myf(b)], edgecolor='k', facecolor='lightgray',linewidth=1)\n",
"\n",
"plt.ylim(bottom=0.0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here we define all three rules"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"def leftpoint(f, a, b, npts):\n",
" '''\n",
" f: Any Python function\n",
" a: Lower integral bound\n",
" b: Upper integral bound\n",
" npts: Number of quadrature points\n",
" \n",
" Returns the integral of f(x) based on the leftpoint rule\n",
" ''' \n",
" x = np.linspace(a,b,npts)\n",
" sum = 0.0\n",
" for i in range(0,len(x)-1):\n",
" a = x[i]\n",
" b = x[i+1]\n",
" sum += (b-a)*(f(a) )\n",
" return sum\n",
"\n",
"def rightpoint(f, a, b, npts):\n",
" '''\n",
" f: Any Python function\n",
" a: Lower integral bound\n",
" b: Upper integral bound\n",
" npts: Number of quadrature points\n",
" \n",
" Returns the integral of f(x) based on the rightpoint rule\n",
" '''\n",
" x = np.linspace(a,b,npts)\n",
" sum = 0.0\n",
" for i in range(0,len(x)-1):\n",
" a = x[i]\n",
" b = x[i+1]\n",
" sum += (b-a)*(f(b) )\n",
" return sum\n",
"\n",
"\n",
"def midpoint(f, a, b, npts):\n",
" '''\n",
" f: Any Python function\n",
" a: Lower integral bound\n",
" b: Upper integral bound\n",
" npts: Number of quadrature points\n",
" \n",
" Returns the integral of f(x) based on the midpoint rule\n",
" '''\n",
" x = np.linspace(a,b,npts)\n",
" sum = 0.0\n",
" for i in range(0,len(x)-1):\n",
" a = x[i]\n",
" b = x[i+1]\n",
" sum += (b-a)*(f( (a+b)/2 ) )\n",
" return sum"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Example"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Use the left, right, and midpoints rules to compute the integral \n",
"\\begin{equation}\n",
"\\int_0^3 e^{-x^2} \\text{d} x\n",
"\\end{equation}\n",
"Compute the absolute relative approximate error for a set of integration points ranging from 2 to 200. Plot the error in all cases."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"numpts = np.arange(2,50)\n",
"a = 0.0\n",
"b = 3.0\n",
"\n",
"lpvals = []\n",
"rpvals = []\n",
"mpvals = []\n",
"for npts in numpts:\n",
" v = leftpoint(myf,a,b,npts) \n",
" lpvals.append(v)\n",
" v = rightpoint(myf,a,b,npts)\n",
" rpvals.append(v)\n",
" v = midpoint(myf,a,b,npts)\n",
" mpvals.append(v)\n",
" \n",
"plt.plot(numpts,lpvals, numpts,rpvals, numpts, mpvals)\n",
"plt.grid()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now let's look at how the absolute approximate relative error behaves"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"lpvals = np.array(lpvals)\n",
"e1 =abs( (lpvals[1:] - lpvals[:-1])/lpvals[1:] )\n",
"plt.loglog(numpts[1:],e1, label=\"Leftpoint rule\")\n",
"\n",
"rpvals = np.array(rpvals)\n",
"e2 =abs( (rpvals[1:] - rpvals[:-1])/rpvals[1:])\n",
"plt.loglog(numpts[1:],e2, label=\"Rightpoint rule\" )\n",
"\n",
"mpvals = np.array(mpvals)\n",
"e3 =abs( (mpvals[1:] - mpvals[:-1])/mpvals[1:])\n",
"plt.loglog(numpts[1:],e3, label=\"Midpoint rule\" )\n",
"\n",
"plt.grid()\n",
"plt.legend()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Trapezoid Rule"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The trapezoidal rule fits a straight line between each two points. It is given by the formula\n",
"\\begin{equation}\n",
"\\int_a^b f(x) \\text{d} x \\approx \\frac{b-a}{2}[f(a) + f(b)]\n",
"\\end{equation}"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0.0, 0.5923990627251441)"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"def myf(x):\n",
" return 1.0/np.sqrt(np.pi) * np.exp(-x**2)\n",
"\n",
"xalot = np.linspace(0,2,200)\n",
"\n",
"npts = 5\n",
"\n",
"x = np.linspace(0,2, npts)\n",
"\n",
"plt.plot(xalot, myf(xalot))\n",
"plt.title('Trapezoidal Rule')\n",
"# ax1.fill_between([a,b],[myf(a),myf(a)], edgecolor='k', facecolor='lightgray',linewidth=1)\n",
"# plot left Riemann Sum\n",
"for i in range(0, len(x)-1):\n",
" a = x[i]\n",
" b = x[i+1]\n",
" plt.fill_between([a,b],[myf(a),myf(b)], edgecolor='k', facecolor='lightgray',linewidth=1)\n",
"\n",
"plt.ylim(bottom=0.0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [],
"source": [
"def traprule(x,y):\n",
" '''\n",
" Integrates a discrete set of data using the Trapezoidal rule.\n",
" \n",
" x: Values of the independent variable\n",
" y: Values of the dependent variable\n",
" \n",
" Returns the integral of f(x) based on the midpoint rule\n",
" '''\n",
"\n",
" sum = 0.0\n",
" for i in range(0,len(x)-1):\n",
" sum += 0.5*(x[i+1]-x[i])*(y[i+1] + y[i])\n",
" return sum\n",
"\n",
"\n",
"def traprulef(f, a, b, npts):\n",
" '''\n",
" Computes the integral of f(x) using the Trapezoidal rule. \n",
" \n",
" f: Any Python function\n",
" a: Lower integral bound\n",
" b: Upper integral bound\n",
" npts: Number of quadrature points\n",
" \n",
" Returns the integral of f(x) based on the midpoint rule\n",
" ''' \n",
" x = np.linspace(a,b,npts)\n",
" sum = 0.0\n",
" for i in range(0,len(x)-1):\n",
" a = x[i]\n",
" b = x[i+1]\n",
" sum += 0.5*(b-a)*(f(a) + f(b))\n",
" return sum"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Example\n",
"Use the Trapezoidal rule to compute the integral \n",
"\\begin{equation}\n",
"\\int_0^3 e^{-x^2} \\text{d} x\n",
"\\end{equation}\n",
"Compute the absolute relative approximate error for a set of integration points ranging from 2 to 200. Plot the error."
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"numpts = np.arange(2,30)\n",
"a = 0.0\n",
"b = 3.0\n",
"\n",
"trapvals = []\n",
"for npts in numpts:\n",
" v = traprulef(myf,a,b,npts) \n",
" trapvals.append(v)\n",
" \n",
"plt.plot(numpts,trapvals)\n",
"plt.grid()"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"trapvals = np.array(trapvals)\n",
"e1 =abs( (trapvals[1:] - trapvals[:-1])/trapvals[1:] )\n",
"plt.loglog(numpts[1:],e1, label=\"Trapezoidal rule\")\n",
"plt.grid()"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"nbpresent": {
"id": "17b00c90-ab55-4e22-a85f-2f09e710ca04"
}
},
"outputs": [
{
"data": {
"text/html": [
"CSS style adapted from https://github.com/barbagroup/CFDPython. Copyright (c) Barba group\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from IPython.core.display import HTML\n",
"def css_styling():\n",
" styles = open(\"../../styles/custom.css\", \"r\").read()\n",
" return HTML(styles)\n",
"css_styling()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"anaconda-cloud": {},
"hide_input": false,
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.7"
}
},
"nbformat": 4,
"nbformat_minor": 2
}