## Numerical Integration: Introduction

### Introduction

Integrals arise naturally in the fields of engineering to describe quantities that are functions of infinitesimal data. For example, if the velocity of a moving object for a particular period of time is known, then the change of position of the moving object can be estimated by summing the velocity at discrete time points multiplied by the corresponding time intervals between the discrete points. The accuracy of the change in position can be increased by decreasing the number of discrete time points, i.e., by reducing the time intervals between the discrete points. The formal definition of an integral of a function is the signed area of the region under the curve between the points and . The definite integral is denoted by:

#### Fundamental Theorem of Calculus

The fundamental theorem of calculus links the concepts of integration and differentiation; roughly speaking, they are the converse of each other. The fundamental theorem of calculus can be stated as follows:

##### Part I

Let be continuous, and let be the function defined such that :

then, is uniformly continuous on , is differentiable on , and is the derivative of (also stated as is the antiderivative of ):

In other words, the first part asserts that if is continuous on the interval , then, an antiderivative always exists; is the derivative of a function that can be defined as with .

##### Part II

Let and be two continuous functions such that (i.e., is the derivative of , or is an antiderivative of ), then:

The second part states that if has an antiderivative , then can be used to calculate the definite integral of on the interval .

For visual proofs of the fundamental theorem of calculus, check out this page.

#### Numerical Integration

The fundamental theorem of calculus allows the calculation of using the antiderivative of . However, if the antiderivative is not available, a numerical procedure can be employed to calculate the integral or the area under the curve. Historically, the term “Quadrature” was used to mean calculating areas. For example, the Gauss numerical integration scheme that will be studied later is also called the Gauss quadrature. The basic problem in numerical integration is to calculate an approximate solution to the definite integral of a function .

##### Riemann Integral

One of the early definitions of the integral of a function is the limit:

where , , and is an arbitrary point such that . If is the sum of the areas of vertical rectangles arranged next to each other and whose top sides touch the function at arbitrary points, then, the Riemann integral is the limit of this sum as the maximum width of the vertical rectangles approaches zero. A function is “Riemann integrable” if the limit above (the limit of ) exists as the width of the rectangles gets smaller and smaller. In this course, we are always dealing with continuous functions. For continuous functions, the limit always exists and thus, all continuous functions are “Riemann integrable”.

##### Newton-Cotes Formulas

The Newton-Cotes formulas rely on replacing the function or tabulated data with an interpolating polynomial that is easy to integrate. In this case, the required integral is evaluated as:

where is a polynomial of degree which can be constructed to pass through data points in the interval .

##### Example

Consider the function defined as:

Calculate the exact integral:

Using steps sizes of , fit an interpolating polynomial to the values of the function at the generated points and calculate the integral by integrating the polynomial. Compare the exact with the interpolating polynomial.

##### Solution

The exact integral can be calculated using Mathematica to be:

The following figures show the interpolating polynomial fitted through data points at step sizes of 0.1, 0.2, 0.5, and 1. The exact integral compared with the numerical integration scheme evaluated by integrating the interpolating polynomial is shown under each curve. Naturally, the higher the degree of the polynomial, the closer the numerical value to the exact value.

The resulting interpolating polynomials are given by:

Integrating the above formulas is straightforward and the resulting values are:

The following is the Mathematica code used.

View Mathematica Code

h = {0.1, 0.2, 0.5, 1}; n = Table[1/h[[i]] + 1, {i, 1, 4}]; y = 2 - x^2 + 0.1*Cos[2 Pi*x/0.7]; Data = Table[{h[[i]]*(j - 1), (y /. x -> h[[i]]*(j - 1.0))}, {i, 1, 4}, {j, 1, n[[i]]}]; Data // MatrixForm; p = Table[InterpolatingPolynomial[Data[[i]], x], {i, 1, 4}]; a1 = Plot[{y, p[[1]]}, {x, 0, 1}, Epilog -> {PointSize[Large], Point[Data[[1]]]}, AxesLabel -> {"x", "y"}, BaseStyle -> Bold, ImageSize -> Medium, PlotRange -> All, Filling -> {2 -> Axis}, PlotLegends -> {"y actual", "y Interpolating Polynomial"}]; a2 = Plot[{y, p[[2]]}, {x, 0, 1}, Epilog -> {PointSize[Large], Point[Data[[2]]]}, AxesLabel -> {"x", "y"}, BaseStyle -> Bold, ImageSize -> Medium, PlotRange -> All, Filling -> {2 -> Axis}, PlotLegends -> {"y actual", "y Interpolating Polynomial"}]; a3 = Plot[{y, p[[3]]}, {x, 0, 1}, Epilog -> {PointSize[Large], Point[Data[[3]]]}, AxesLabel -> {"x", "y"}, BaseStyle -> Bold, ImageSize -> Medium, PlotRange -> All, Filling -> {2 -> Axis}, PlotLegends -> {"y actual", "y Interpolating Polynomial"}]; a4 = Plot[{y, p[[4]]}, {x, 0, 1}, Epilog -> {PointSize[Large], Point[Data[[4]]]}, AxesLabel -> {"x", "y"}, BaseStyle -> Bold, ImageSize -> Medium, PlotRange -> All, Filling -> {2 -> Axis}, PlotLegends -> {"y actual", "y Interpolating Polynomial"}]; Iexact = Integrate[y, {x, 0, 1}]; Inumerical = Table[Integrate[p[[i]], {x, 0, 1}], {i, 1, 4}]; atext = Table[Grid[{{"I exact = ", Iexact, ", I numerical = ", Inumerical[[i]]}}], {i, 1, 4}]; Grid[{{a1, a2}, {atext[[1]], atext[[2]]}}] Grid[{{a3, a4}, {atext[[3]], atext[[4]]}}]

import numpy as np import sympy as sp import matplotlib.pyplot as plt h = [0.1, 0.2, 0.5, 1] n = [1/h[i] + 1 for i in range(len(h))] x = sp.symbols('x') y = 2 - x**2 + 0.1*sp.cos(2*np.pi*x/0.7) Data = [[(h[i]*j, (y.subs(x, h[i]*j))) for j in range(int(n[i]))] for i in range(len(h))] p = [sp.interpolate(Data[i], x) for i in range(len(h))] x_val = np.arange(0,1,0.01) Iexact = sp.integrate(y, (x, 0, 1)) plt.xlabel("I exact = "+str(sp.N(Iexact,6)) + " , I numerical = "+str(sp.N(sp.integrate(p[0], (x, 0, 1)),6))) plt.scatter([point[0] for point in Data[0]],[point[1] for point in Data[0]], c='k') plt.plot(x_val, [y.subs(x,i) for i in x_val], label="y actual") plt.plot(x_val, [p[0].subs(x,i) for i in x_val], label="y Interpolating Polynomial") plt.fill_between(x_val, np.float_([p[0].subs(x,i) for i in x_val]), 0, alpha=0.20) plt.legend(); plt.grid(); plt.show() plt.xlabel("I exact = "+str(sp.N(Iexact,6)) + " , I numerical = "+str(sp.N(sp.integrate(p[1], (x, 0, 1)),6))) plt.scatter([point[0] for point in Data[1]],[point[1] for point in Data[1]], c='k') plt.plot(x_val, [y.subs(x,i) for i in x_val], label="y actual") plt.plot(x_val, [p[1].subs(x,i) for i in x_val], label="y Interpolating Polynomial") plt.fill_between(x_val, np.float_([p[1].subs(x,i) for i in x_val]), 0, alpha=0.20) plt.legend(); plt.grid(); plt.show() plt.xlabel("I exact = "+str(sp.N(Iexact,6)) + " , I numerical = "+str(sp.N(sp.integrate(p[2], (x, 0, 1)),6))) plt.scatter([point[0] for point in Data[2]],[point[1] for point in Data[2]], c='k') plt.plot(x_val, [y.subs(x,i) for i in x_val], label="y actual") plt.plot(x_val, [p[2].subs(x,i) for i in x_val], label="y Interpolating Polynomial") plt.fill_between(x_val, np.float_([p[2].subs(x,i) for i in x_val]), 0, alpha=0.20) plt.legend(); plt.grid(); plt.show() plt.xlabel("I exact = "+str(sp.N(Iexact,6)) + " , I numerical = "+str(sp.N(sp.integrate(p[3], (x, 0, 1)),6))) plt.scatter([point[0] for point in Data[3]],[point[1] for point in Data[3]], c='k') plt.plot(x_val, [y.subs(x,i) for i in x_val], label="y actual") plt.fill_between(x_val, np.float_([p[3].subs(x,i) for i in x_val]), 0, alpha=0.20) plt.plot(x_val, [p[3].subs(x,i) for i in x_val], label="y Interpolating Polynomial") plt.legend(); plt.grid(); plt.show()

The following links provide the MATLAB codes for implementing the integration using the interpolating polynomials

However, as described in the interpolating polynomial section, the larger the degree of the polynomial, the more susceptible it is to oscillations. As an example, the previous problem is repeated with the Runge function defined on the interval from -1 to 1. The figure below shows that using many points leads to large oscillations in the interpolating polynomial which renders the numerical integration largely inaccurate. When fewer points are used, the interpolating polynomial is still highly inaccurate. Therefore, the Newton-Cotes formulas can be applied by simply subdividing the interval into smaller subintervals and applying the Newton-Cotes formulas to each subinterval and adding the results. This process is called: “Composite rules”. The rectangle method, trapezoidal rule, and Simpson’s rules presented in the following sections are specific examples of the composite Newton-Cotes formulas.

h = {0.1, 0.2, 0.5, 1}; n = Table[2/h[[i]] + 1, {i, 1, 4}]; y = 1.0/(1 + 25 x^2); Data = Table[{-1 + h[[i]]*(j - 1), (y /. x -> -1 + h[[i]]*(j - 1.0))}, {i, 1, 4}, {j, 1, n[[i]]}]; Data // MatrixForm; p = Table[InterpolatingPolynomial[Data[[i]], x], {i, 1, 4}]; a1 = Plot[{y, p[[1]]}, {x, -1, 1}, Epilog -> {PointSize[Large], Point[Data[[1]]]}, AxesLabel -> {"x", "y"}, BaseStyle -> Bold, ImageSize -> Medium, PlotRange -> All, Filling -> {2 -> Axis}, PlotLegends -> {"y actual", "y Interpolating Polynomial"}]; a2 = Plot[{y, p[[2]]}, {x, -1, 1}, Epilog -> {PointSize[Large], Point[Data[[2]]]}, AxesLabel -> {"x", "y"}, BaseStyle -> Bold, ImageSize -> Medium, PlotRange -> All, Filling -> {2 -> Axis}, PlotLegends -> {"y actual", "y Interpolating Polynomial"}]; a3 = Plot[{y, p[[3]]}, {x, -1, 1}, Epilog -> {PointSize[Large], Point[Data[[3]]]}, AxesLabel -> {"x", "y"}, BaseStyle -> Bold, ImageSize -> Medium, PlotRange -> All, Filling -> {2 -> Axis}, PlotLegends -> {"y actual", "y Interpolating Polynomial"}]; a4 = Plot[{y, p[[4]]}, {x, -1, 1}, Epilog -> {PointSize[Large], Point[Data[[4]]]}, AxesLabel -> {"x", "y"}, BaseStyle -> Bold, ImageSize -> Medium, PlotRange -> All, Filling -> {2 -> Axis}, PlotLegends -> {"y actual", "y Interpolating Polynomial"}]; Iexact = Integrate[y, {x, -1, 1}]; Inumerical = Table[Integrate[p[[i]], {x, -1, 1}], {i, 1, 4}]; atext = Table[ Grid[{{"I exact = ", Iexact, ", I numerical = ", Inumerical[[i]]}}], {i, 1, 4}]; Grid[{{a1, a2}, {atext[[1]], atext[[2]]}}] Grid[{{a3, a4}, {atext[[3]], atext[[4]]}}]

import numpy as np import sympy as sp import matplotlib.pyplot as plt h = [0.1, 0.2, 0.5, 1] n = [2/h[i] + 1 for i in range(len(h))] x = sp.symbols('x') y = 1.0/(1 + 25*x**2) Data = [[(-1 + h[i]*j, (y.subs(x, -1 + h[i]*j))) for j in range(int(n[i]))] for i in range(len(h))] p = [sp.interpolate(Data[i], x) for i in range(len(h))] x_val = np.arange(-1,1,0.01) Iexact = sp.integrate(y, (x, -sp.oo, sp.oo)) plt.xlabel("I exact = "+str(sp.N(Iexact,6)) + " , I numerical = "+str(sp.N(sp.integrate(p[0], (x, -1, 1)),6))) plt.scatter([point[0] for point in Data[0]],[point[1] for point in Data[0]], c='k') plt.plot(x_val, [y.subs(x,i) for i in x_val], label="y actual") plt.plot(x_val, [p[0].subs(x,i) for i in x_val], label="y Interpolating Polynomial") plt.fill_between(x_val, np.float_([p[0].subs(x,i) for i in x_val]), 0, alpha=0.20) plt.legend(); plt.grid(); plt.show() plt.xlabel("I exact = "+str(sp.N(Iexact,6)) + " , I numerical = "+str(sp.N(sp.integrate(p[1], (x, -1, 1)),6))) plt.scatter([point[0] for point in Data[1]],[point[1] for point in Data[1]], c='k') plt.plot(x_val, [y.subs(x,i) for i in x_val], label="y actual") plt.plot(x_val, [p[1].subs(x,i) for i in x_val], label="y Interpolating Polynomial") plt.fill_between(x_val, np.float_([p[1].subs(x,i) for i in x_val]), 0, alpha=0.20) plt.legend(); plt.grid(); plt.show() plt.xlabel("I exact = "+str(sp.N(Iexact,6)) + " , I numerical = "+str(sp.N(sp.integrate(p[2], (x, -1, 1)),6))) plt.scatter([point[0] for point in Data[2]],[point[1] for point in Data[2]], c='k') plt.plot(x_val, [y.subs(x,i) for i in x_val], label="y actual") plt.plot(x_val, [p[2].subs(x,i) for i in x_val], label="y Interpolating Polynomial") plt.fill_between(x_val, np.float_([p[2].subs(x,i) for i in x_val]), 0, alpha=0.20) plt.legend(); plt.grid(); plt.show() plt.xlabel("I exact = "+str(sp.N(Iexact,6)) + " , I numerical = "+str(sp.N(sp.integrate(p[3], (x, -1, 1)),6))) plt.scatter([point[0] for point in Data[3]],[point[1] for point in Data[3]], c='k') plt.plot(x_val, [y.subs(x,i) for i in x_val], label="y actual") plt.fill_between(x_val, np.float_([p[3].subs(x,i) for i in x_val]), 0, alpha=0.20) plt.plot(x_val, [p[3].subs(x,i) for i in x_val], label="y Interpolating Polynomial") plt.legend(); plt.grid(); plt.show()