What is Numerical Integration?

In a mathematical investigation, numerical integration comprises a wide group of calculations for computing the mathematical estimation of definite integral, and likewise, the term is moreover in some cases used to depict the numerical solution of differential equations. The most important aspect of this theory is error analysis.

Important Aspects of Numerical Integration

The term was coined in 1656, by the mathematician John Wallis in his series ‘Arithmetical Infinitorum’.

The term is abbreviated as quadrature. The set of techniques to execute the process of integration are grouped under Newton-Cotes formulas. The aim is to evaluate the integrand using midpoint, trapezoidal, and Simpson’s rule.

When do we use Numerical Integration?

• Sometimes, we may find that the integrand is either regularly or irregularly spaced or some data values are missing when we are working on experimental data.
• The upper or lower limit of integration is infinite. There are some points where the function is discontinuous.

How Numerical Integration is different from Analytical Integration?

• The analytic integration is computed below,

The motivation behind introducing the concept of mathematics is some integrals like ${\displaystyle \underset{8}{\overset{12}{\int }}{e}^{-x^3}dx}$.

• Analytical integration always gives an exact answer whereas, numerical integration is an approximation of definite integral.
• This only gives the approximate area under the curve.
• The objective is to find the approximate area between the graph of continuous function and the x-axis.

Newton-Cotes Formula

We utilize the concept of interpolation in order to integrate the function and it is approximated by $n^{\text{th}}$ degree polynomial.

Consider a function $f\left(x\right)$ that is approximated by a polynomial over a closed interval $\left[a,b\right]$.

${\displaystyle \underset{a}{\overset{b}{\int }}f\left(x\right)dx\approx }{\displaystyle \underset{a}{\overset{b}{\int }}\left(a_0+a_{1}x+\cdots +a_n{x}^n\right)dx}$

On further simplification we get,

$\begin{array}{c}{\displaystyle \underset{a}{\overset{b}{\int }}f\left(x\right)dx}={\left[a_{0}x+\frac{a_1{x}^2}{2}+\cdots +\frac{a_n{x}^{n+1}}{n+1}\right]}_a^b\\ =a_0\left(b-a\right)+\frac{a_1\left(b^2-a^2\right)}{2}+\cdots +\frac{a_n\left(b^{n+1}-a^{n+1}\right)}{n+1}\end{array}$

Euler’s Method

The Euler method (also known as the forward Euler method) is a first-order numerical technique for solving ODEs with a given initial value. It's the most fundamental explicit approach for numerically integrating ordinary differential equations.

We use Euler's Method to produce a numerical solution to a problem with an initial value of the form:

$y\text{'}=f\left(x,y\right)$

$y\left(x_0\right)=y_0$

We decide on what interval we want to find the solution in, starting with the initial condition. This interval is divided into h-sized subdivisions. Then, using the iterative formulas, we produce the rest of the solution, beginning with the initial condition:

$x_{n+1}=x_n+h$

$y_{n+1}=y_n+hf\left(x_n,y_n\right)$

To determine the points' coordinates in our numerical solution when we hit the right end of the desired interval, we stop the process.

Techniques of Numerical Integration

Midpoint Rule

Let a function $g\left(x\right)$ be defined on the closed and continuous interval $\left[a,b\right]$. We partition $\left[a,b\right]$in form of rectangles into n subintervals where the length of each subinterval is denoted by $\Delta x$and mathematically this length is expressed as $\Delta x=\frac{b-a}{n}$.

For an arbitrary $i^{\text{th}}$ interval if $m_i$ is the midpoint, then midpoint rule is expressed as,

$M_n={\displaystyle \sum _{i=1}^{n}g\left(m_i\right)}\Delta x$

As the number of rectangles becomes large in number, that is $n\to \infty$, the definite integral is expressed as:

$\underset{n\to \infty }{\lim }{M}_n={\displaystyle \underset{a}{\overset{b}{\int }}g\left(x\right)dx}$

The definite error bound for the midpoint rule is expressed as:

$M_n\le \frac{M{\left(b-a\right)}^3}{24n^2}$

Since we don’t use the endpoints of $\left[a,b\right]$in computing the value, this method is known as an open method.

Trapezoidal Rule

Let a function $g\left(x\right)$ be defined on the closed and continuous interval $\left[a,b\right]$. We partition $\left[a,b\right]$in form of trapezoid into n subintervals where the length of each subinterval is denoted by $\Delta x$and mathematically this length is expressed as $\Delta x=\frac{b-a}{n}$.

We can label the points on the x-axis as $P=\left\{x_0,x_1,x_2,\dots ,x_n\right\}$.

The area of a trapezoid is half times the sum of parallel side times the distance between them.

We can express the area of the first trapezoid as,

$A_1=\frac{\Delta x}{2}\left[g\left(x_0\right)+g\left(x_1\right)\right]$

Similarly, the area of the second trapezoid is expressed as,

$A_2=\frac{\Delta x}{2}\left[g\left(x_1\right)+g\left(x_2\right)\right]$

And continuing the pattern, the area of the $n^{\text{th}}$ trapezoid is expressed as,

$A_n=\frac{\Delta x}{2}\left[g\left(x_{n-1}\right)+g\left(x_n\right)\right]$

Summing up the area of all the trapezoids, the trapezoidal rule approximation is expressed as,

$T_n=\frac{\Delta x}{2}\left[g\left(x_0\right)+2g\left(x_1\right)+2g\left(x_2\right)+\cdots +2g\left(x_{n-1}\right)+g\left(x_n\right)\right]$

As the number of trapezoids becomes large in number, that is $n\to \infty$, the definite integral is expressed as,

$\underset{n\to \infty }{\lim }{T}_n={\displaystyle \underset{a}{\overset{b}{\int }}g\left(x\right)dx}$

The expression that gives an error bound for the trapezoidal rule is provided below.

$T_n\le \frac{M{\left(b-a\right)}^3}{12n^2}$

Here M is the maximum value of the second derivative of the function $g\left(x\right)$ over $\left[a,b\right]$, that is $M=\max \left|g\text{'}\text{'}\left(x\right)\right|$.

Let’s practice an example based on this methodology.

Problem 1:

The objective is to compute the area under the curve of the function $g\left(x\right)=\frac{1}{2+x}$ over $\left[0,1\right]$ and partitioning $\left[0,1\right]$into four subintervals using trapezoidal rule approximation.

Answer: First, we partition $\left[0,1\right]$into four subintervals of equal length.

The length of each interval is expressed as $\Delta x=\frac{b-a}{n}$.

Substituting b as 1, a as 0, and n as 4 in the above expression,

$\begin{array}{c}\Delta x=\frac{1-0}{4}\\ =0.25\end{array}$

Therefore, the partition is $P=\left\{0,0.25,0.5,0.75,1\right\}$.

The sub-intervals are: $\left[0,0.25\right],\left[0.25,0.5\right],\left[0.5,0.75\right]\text{ and }\left[0.75,1\right]$.

Let $y=\frac{1}{2+x}$. Construct a table of values by estimating the value of the function at the endpoints of the sub-intervals of $\left[0,1\right]$.

Applying the trapezoidal rule we have,

${\displaystyle \int g\left(x\right)dx}=\frac{\Delta x}{2}\left[y_0+2\left(y_1+y_2+y_3\right)+y_4\right]$

$=\frac{0.25}{2}\left[0.5+2\left(0.4444+0.4+0.3636\right)+0.3333\right]$

$=\frac{0.25}{2}\left[0.5+2\left(1.2081\right)+0.3333\right]$

$=0.4062$

Thus, the approximate area under the curve of the function $g\left(x\right)=\frac{1}{2+x}$ using trapezoidal rule over $\left[0,1\right]$ is 0.4062 square units.

Simpson’s One-Third Rule

Let a function $g\left(x\right)$ be defined on the closed and continuous interval $\left[a,b\right]$. We partition this interval in form of piecewise quadratic functions into n subintervals where the length of each subinterval is denoted by $\Delta x$and mathematically this length is expressed as $\Delta x=\frac{b-a}{n}$.

Summing up the area of all the quadratic functions, Simpson’s one-third rule approximation is expressed as,

$S_n=\frac{\Delta x}{3}\left[f\left(x_0\right)+4f\left(x_1\right)+2f\left(x_2\right)+4f\left(x_3\right)+\cdots +2f\left(x_{n-2}\right)+4f\left(x_{n-1}\right)+f\left(x_n\right)\right]$

Similarly, we can express Simpson’s three-eighth rule approximation as,

$S_n=\frac{3\Delta x}{8}\left[f\left(x_0\right)+3f\left(x_1\right)+3f\left(x_2\right)+2f\left(x_3\right)+\cdots +3f\left(x_{n-2}\right)+3f\left(x_{n-1}\right)+f\left(x_n\right)\right]$

As the number of quadratic functions becomes large in number, that is $n\to \infty$, the definite integral is expressed as,

$\underset{n\to \infty }{\lim }{S}_n={\displaystyle \underset{a}{\overset{b}{\int }}g\left(x\right)dx}$

The expression that gives an error bound for Simpson’s rule is provided below.

$S_n\le \frac{M{\left(b-a\right)}^5}{180n^4}$

Here M is the maximum value of the fourth derivative of the function $g\left(x\right)$ over $\left[a,b\right]$, that is $M=\max \left|g^{\left(4\right)}\left(x\right)\right|$.

Since we use the endpoints of the interval in computing the value, this method is known as a closed method.

Problem 2:

The objective is to compute the area under the curve of the function $g\left(x\right)=\frac{1}{x^2+2x+3}$ over $\left[0,1\right]$ and partitioning the interval into four subintervals using Simpson’s one-third rule approximation.

Answer: First we partition the interval $\left[0,1\right]$ into four subintervals of equal length.

The length of each interval is expressed as $\Delta x=\frac{b-a}{n}$.

Substituting b as 1, a as 0, and n as 4 in the above expression,

$\begin{array}{c}\Delta x=\frac{1-0}{4}\\ =0.25\end{array}$

Therefore, the $\left[0,1\right]$is partitioned as $P=\left\{0,0.25,0.5,0.75,1\right\}$.

The sub-intervals are: $\left[0,0.25\right],\left[0.25,0.5\right],\left[0.5,0.75\right]\text{ and }\left[0.75,1\right]$.

Let $g\left(x\right)=\frac{1}{x^2+2x+3}$. Construct a table of values by estimating the value of the function at the endpoints of the sub-intervals of $\left[0,1\right]$.

Applying Simpson’s rule we have,

${\displaystyle \int g\left(x\right)dx}=\frac{\Delta x}{3}\left[y_0+2\left(y_2\right)+4\left(y_1+y_3\right)+y_4\right]$

$=\frac{0.25}{2}\left[0.3333+2\left(0.2353\right)+4\left(0.2807+0.1975\right)+0.1667\right]$

$=0.2403$

Thus, the approximate area under the curve of the function $g\left(x\right)=\frac{1}{x^2+2x+3}$ over $\left[0,1\right]$ is 0.2403square units.

Problem 3:

The objective is to compute the area under the curve of the function $g\left(x\right)=3x^3-5x+1$ over the interval $\left[2,4\right]$ and partitioning it into four subintervals using Simpson’s three-eighth rule approximation.

Answer: First, we partition $\left[2,4\right]$into four sub-intervals of equal length.

The length of each interval is expressed as $\Delta x=\frac{b-a}{n}$.

Substituting b as 4, a as 2, and n as 4 in the above expression,

$\begin{array}{c}\Delta x=\frac{4-2}{4}\\ =0.5\end{array}$

Therefore, the interval $\left[2,4\right]$is partitioned as $P=\left\{2,2.5,3,3.5,4\right\}$.

The sub-intervals are: $\left[2,2.5\right],\left[2.5,3\right],\left[3,3.5\right]\text{ and }\left[3.5,4\right]$.

Let $g\left(x\right)=3x^3-5x+1$. Construct a table of values by estimating the value of the function at the endpoints of the sub-intervals of $\left[2,4\right]$.

Applying the Simpson’s three-eighth rule we have,

${\displaystyle \int g\left(x\right)dx}=\frac{3\Delta x}{8}\left[y_0+2\left(y_3\right)+3\left(y_1+y_2\right)+y_4\right]$

$=\frac{3\left(0.5\right)}{8}\left[15+2\left(112.125\right)+3\left(35.375+67\right)+173\right]$

$=134.8828$

Thus, the approximate area under the curve of the function $g\left(x\right)=3x^3-5x+1$ over $\left[2,4\right]$ is 134.8828 square units.

Formulas

• Newton-Cotes formula for the approximation of a function by a polynomial:

${\displaystyle \underset{a}{\overset{b}{\int }}f\left(x\right)dx\approx }{\displaystyle \underset{a}{\overset{b}{\int }}\left(a_0+a_{1}x+\cdots +a_n{x}^n\right)dx}$

$T_n=\frac{\Delta x}{2}\left[g\left(x_0\right)+2g\left(x_1\right)+2g\left(x_2\right)+\cdots +2g\left(x_{n-1}\right)+g\left(x_n\right)\right]$, where $\underset{n\to \infty }{\lim }{T}_n={\displaystyle \underset{a}{\overset{b}{\int }}g\left(x\right)dx}$.

$S_n=\frac{\Delta x}{3}\left[f\left(x_0\right)+4f\left(x_1\right)+2f\left(x_2\right)+4f\left(x_3\right)+\cdots +2f\left(x_{n-2}\right)+4f\left(x_{n-1}\right)+f\left(x_n\right)\right]$, where $\underset{n\to \infty }{\lim }{S}_n={\displaystyle \underset{a}{\overset{b}{\int }}g\left(x\right)dx}$.

Context and Applications

The demand for the techniques of this field is very high because of the high-tech software available to compute them. We can convert the real-life problem into mathematical expression and apply the techniques of numerical integration for better results to the problems.

• One of the first applications of integration was in determining the volume of wine casks (with curved surfaces). Other applications of integration include finding areas under curved surfaces, finding centers of mass, displacement and velocity, and fluid flow.

• In Geology, to estimate the volume of hydrocarbon reservoirs that are in form of close structural traps the techniques of this field of mathematics are used.

• The main goal of meteorological stations is to create a framework that allows for effective climate change monitoring. The mathematical investigation using these techniques is of utmost importance.

Practice Problem

Question: Compute the approximate value at $x=0.5$for the differential equation $y\text{'}=5x+2y$ with step size as 0.1 and initial condition $y\left(0\right)=3$ using Euler method approximation.

Solution: The objective is to estimate the value of $y\left(0.5\right)$.

We are provided that $f\left(x,y\right)=5x+2y$ and $h=0.1$. Applying the Euler method:

The first iteration is,

$y_1=y_0+hf\left(x_0,y_0\right)$

$=3+0.1f\left(0,3\right)$

$=3+0.1\left(6\right)$

$=3.6$

Second iteration:

$y_2=y_1+hf\left(x_1,y_1\right)$

$=3.6+0.1f\left(0.1,3.6\right)$

$=3.6+0.1\left(7.7\right)$

$=4.37$

Third iteration:

$y_3=y_2+hf\left(x_2,y_2\right)$

$=4.37+0.1f\left(0.2,4.37\right)$

$=4.37+0.1\left(9.74\right)$

$=5.34$

Fourth iteration:

$y_4=y_3+hf\left(x_3,y_3\right)$

$=5.34+0.1f\left(0.3,5.34\right)$

$=5.34+0.1\left(12.19\right)$

$=6.56$

Fifth iteration:

$y_5=y_4+hf\left(x_4,y_4\right)$

$=6.56+0.1f\left(0.4,6.56\right)$

$=6.56+0.1\left(15.13\right)$

Thus, the value at $x=0.5$ is 8.08 that is $y\left(0.5\right)=8.08$.