## Numerical Integration: Gauss Quadrature

### Gauss Quadrature

The Gauss integration scheme is a very efficient method to perform numerical integration over intervals. In fact, if the function to be integrated is a polynomial of an appropriate degree, then the Gauss integration scheme produces exact results. The Gauss integration scheme has been implemented in almost every finite element analysis software due to its simplicity and computational efficiency. This section outlines the basic principles behind the Gauss integration scheme. For its application to the finite element analysis method, please visit this section.

#### Introduction

Gauss quadrature aims to find the “least” number of fixed points to approximate the integral of a function such that:

where and . Also, is called an integration point and is called the associated weight. The number of integration points and the associated weights are chosen according to the complexity of the function to be integrated. Since a general polynomial of degree has coefficients, it is possible to find a Gauss integration scheme with number of integration points and number of associated weights to exactly integrate that polynomial function on the interval . The following figure illustrates the concept of using the Gauss integration points to calculate the area under the curve for polynomials. In the following sections, the required number of integration points for particular polynomials are presented.

#### Affine Functions (First-Degree Polynomials): One Integration Point

For a first-degree polynomial, , therefore, , so, 1 integration point is sufficient as will be shown here. Consider with , then:

So, for functions that are very close to being affine, a numerical integration scheme with 1 integration point that is with an associated weight of 2 can be employed. In other words, a one-point numerical integration scheme has the form:

#### Third-Degree Polynomials: Two Integration Points

For a third-degree polynomial, , therefore, , so, 2 integration points are sufficient as will be shown here. Consider , then:

Assuming 2 integration points in the Gauss integration scheme yields:

For the Gauss integration scheme to yield accurate results, the right-hand sides of the two equations above need to be equal for any choice of a cubic function. Therefore the multipliers of , , , and should be the same. So, four equations in four unknowns can be written as follows:

The solution to the above four equations yields , , . The Mathematica code below forms the above equations and solves them:

View Mathematica Codef = a0 + a1*x + a2*x^2 + a3*x^3; I1 = Integrate[f, {x, -1, 1}] I2 = w1*(f /. x -> x1) + w2*(f /. x -> x2) Eq1 = Coefficient[I1 - I2, a0] Eq2 = Coefficient[I1 - I2, a1] Eq3 = Coefficient[I1 - I2, a2] Eq4 = Coefficient[I1 - I2, a3] Solve[{Eq1 == 0, Eq2 == 0, Eq3 == 0, Eq4 == 0}, {x1, x2, w1, w2}]

import sympy as sp sp.init_printing(use_latex=True) x, x1, x2 = sp.symbols('x x1 x2') a0, a1, a2, a3 = sp.symbols('a0 a1 a2 a3') w1, w2 = sp.symbols('w1 w2') f = a0 + a1*x + a2*x**2 + a3*x**3 I1 = sp.integrate(f, (x, -1, 1)) I2 = w1*(f.subs(x, x1)) + w2*(f.subs(x, x2)) display("I1: ",I1) display("I2: ",I2) Eq1 = sp.expand(I1 - I2).coeff(a0) Eq2 = sp.expand(I1 - I2).coeff(a1) Eq3 = sp.expand(I1 - I2).coeff(a2) Eq4 = sp.expand(I1 - I2).coeff(a3) display("Eq1: ",Eq1) display("Eq2: ",Eq2) display("Eq3: ",Eq3) display("Eq4: ",Eq4) sol = list(sp.nonlinsolve([Eq1, Eq2, Eq3, Eq4], [x1, x2, w1, w2])) display(sol)

So, for functions that are very close to being cubic, the following numerical integration scheme with 2 integration points can be employed:

#### Fifth-Degree Polynomials: Three Integration Points

For a fifth-degree polynomial, , therefore, , so, 3 integration points are sufficient to exactly integrate a fifth-degree polynomial. Consider , then:

Assuming 3 integration points in the Gauss integration scheme yields:

For the Gauss integration scheme to yield accurate results, the right-hand sides of the two equations above need to be equal for any choice of a fifth-order polynomial function. Therefore the multipliers of , , , , , and should be the same. So, six equations in six unknowns can be written as follows:

The solution to the above six equations yields , , , , , and . The Mathematica code below forms the above equations and solves them:

View Mathematica Codef = a0 + a1*x + a2*x^2 + a3*x^3 + a4*x^4 + a5*x^5; I1 = Integrate[f, {x, -1, 1}] I2 = w1*(f /. x -> x1) + w2*(f /. x -> x2) + w3*(f /. x -> x3) Eq1 = Coefficient[I1 - I2, a0] Eq2 = Coefficient[I1 - I2, a1] Eq3 = Coefficient[I1 - I2, a2] Eq4 = Coefficient[I1 - I2, a3] Eq5 = Coefficient[I1 - I2, a4] Eq6 = Coefficient[I1 - I2, a5] Solve[{Eq1 == 0, Eq2 == 0, Eq3 == 0, Eq4 == 0, Eq5 == 0, Eq6 == 0}, {x1, x2, x3, w1, w2, w3}]

import sympy as sp sp.init_printing(use_latex=True) x, x1, x2, x3 = sp.symbols('x x1 x2 x3') a0, a1, a2, a3, a4, a5 = sp.symbols('a0 a1 a2 a3 a4 a5') w1, w2, w3 = sp.symbols('w1 w2 w3') f = a0 + a1*x + a2*x**2 + a3*x**3 + a4*x**4 + a5*x**5 I1 = sp.integrate(f, (x, -1, 1)) I2 = w1*(f.subs(x, x1)) + w2*(f.subs(x, x2)) + w3*(f.subs(x, x3)) display("I1: ",I1) display("I2: ",I2) Eq1 = sp.expand(I1 - I2).coeff(a0) Eq2 = sp.expand(I1 - I2).coeff(a1) Eq3 = sp.expand(I1 - I2).coeff(a2) Eq4 = sp.expand(I1 - I2).coeff(a3) Eq5 = sp.expand(I1 - I2).coeff(a4) Eq6 = sp.expand(I1 - I2).coeff(a5) display("Eq1: ",Eq1) display("Eq2: ",Eq2) display("Eq3: ",Eq3) display("Eq4: ",Eq4) display("Eq5: ",Eq5) display("Eq6: ",Eq6) sol = list(sp.nonlinsolve([Eq1, Eq2, Eq3, Eq4, Eq5, Eq6], [x1, x2, x3, w1, w2, w3])) display(sol)

So, for functions that are very close to being fifth-order polynomials, the following numerical integration scheme with 3 integration points can be employed:

#### Seventh-Degree Polynomial: Four Integration Points

For a seventh-degree polynomial, , therefore, , so, 4 integration points are sufficient to exactly integrate a seventh-degree polynomial. Following the procedure outlined above, the integration points along with the associated weight:

So, for functions that are very close to being seventh-order polynomials, the following numerical integration scheme with 4 integration points can be employed:

The following is the Mathematica code used to find the above integration points and corresponding weights.

View Mathematica Codef = a0 + a1*x + a2*x^2 + a3*x^3 + a4*x^4 + a5*x^5 + a6*x^6 + a7*x^7; I1 = Integrate[f, {x, -1, 1}] I2 = w1*(f /. x -> x1) + w2*(f /. x -> x2) + w3*(f /. x -> x3) + w4*(f /. x -> x4) Eq1 = Coefficient[I1 - I2, a0] Eq2 = Coefficient[I1 - I2, a1] Eq3 = Coefficient[I1 - I2, a2] Eq4 = Coefficient[I1 - I2, a3] Eq5 = Coefficient[I1 - I2, a4] Eq6 = Coefficient[I1 - I2, a5] Eq7 = Coefficient[I1 - I2, a6] Eq8 = Coefficient[I1 - I2, a7] a = Solve[{Eq1 == 0, Eq2 == 0, Eq3 == 0, Eq4 == 0, Eq5 == 0, Eq6 == 0, Eq7 == 0, Eq8 == 0}, {x1, x2, x3, w1, w2, w3, x4, w4}] Chop[N[a[[1]]]]

import sympy as sp sp.init_printing(use_latex=True) x, x1, x2, x3, x4 = sp.symbols('x x1 x2 x3 x4') a0, a1, a2, a3, a4, a5, a6, a7 = sp.symbols('a0 a1 a2 a3 a4 a5 a6 a7') w1, w2, w3, w4 = sp.symbols('w1 w2 w3 w4') f = a0 + a1*x + a2*x**2 + a3*x**3 + a4*x**4 + a5*x**5 + a6*x**6 + a7*x**7 I1 = sp.integrate(f, (x, -1, 1)).subs({x3:-x1,x4:-x2,w3:w1,w4:w2}) I2 = (w1*(f.subs(x, x1)) + w2*(f.subs(x, x2)) + w3*(f.subs(x, x3)) + w4*(f.subs(x, x4))).subs({x3:-x1,x4:-x2,w3:w1,w4:w2}) display("I1: ",I1) display("I2: ",I2) Eq1 = sp.expand(I1 - I2).coeff(a0) Eq2 = sp.expand(I1 - I2).coeff(a1) Eq3 = sp.expand(I1 - I2).coeff(a2) Eq4 = sp.expand(I1 - I2).coeff(a3) Eq5 = sp.expand(I1 - I2).coeff(a4) Eq6 = sp.expand(I1 - I2).coeff(a5) Eq7 = sp.expand(I1 - I2).coeff(a6) Eq8 = sp.expand(I1 - I2).coeff(a7) display("Eq1: ",Eq1) display("Eq2: ",Eq2) display("Eq3: ",Eq3) display("Eq4: ",Eq4) display("Eq5: ",Eq5) display("Eq6: ",Eq6) display("Eq7: ",Eq7) display("Eq8: ",Eq8) sol = list(sp.nonlinsolve([Eq1, Eq3, Eq5, Eq7], [x1, x2, w1, w2])) display(sol)

#### Ninth-Degree Polynomial: Five Integration Points

For a ninth-degree polynomial, , therefore, , so, 5 integration points are sufficient to exactly integrate a ninth-degree polynomial. Following the procedure outlined above, the integration points along with the associated weight:

So, for functions that are very close to being ninth-order polynomials, the following numerical integration scheme with 5 integration points can be employed:

The following is the Mathematica code used to find the above integration points and corresponding weights.

View Mathematica Codef = a0 + a1*x + a2*x^2 + a3*x^3 + a4*x^4 + a5*x^5 + a6*x^6 + a7*x^7 + a8*x^8 + a9*x^9; I1 = Integrate[f, {x, -1, 1}] I2 = w1*(f /. x -> x1) + w2*(f /. x -> x2) + w3*(f /. x -> x3) + w4*(f /. x -> x4) + w5*(f /. x -> x5) Eq1 = Coefficient[I1 - I2, a0] Eq2 = Coefficient[I1 - I2, a1] Eq3 = Coefficient[I1 - I2, a2] Eq4 = Coefficient[I1 - I2, a3] Eq5 = Coefficient[I1 - I2, a4] Eq6 = Coefficient[I1 - I2, a5] Eq7 = Coefficient[I1 - I2, a6] Eq8 = Coefficient[I1 - I2, a7] Eq9 = Coefficient[I1 - I2, a8] Eq10 = Coefficient[I1 - I2, a9] a = Solve[{Eq1 == 0, Eq2 == 0, Eq3 == 0, Eq4 == 0, Eq5 == 0, Eq6 == 0, Eq7 == 0, Eq8 == 0, Eq9 == 0, Eq10 == 0}, {x1, x2, x3, w1, w2, w3, x4, w4, x5, w5}] Chop[N[a[[1]]]]

import sympy as sp sp.init_printing(use_latex=True) x, x1, x2, x3, x4, x5 = sp.symbols('x x1 x2 x3 x4 x5') a0, a1, a2, a3, a4, a5, a6, a7, a8, a9 = sp.symbols('a0 a1 a2 a3 a4 a5 a6 a7 a8 a9') w1, w2, w3, w4, w5 = sp.symbols('w1 w2 w3 w4 w5') f = a0 + a1*x + a2*x**2 + a3*x**3 + a4*x**4 + a5*x**5 + a6*x**6 + a7*x**7 + a8*x**8 + a9*x**9 I1 = sp.integrate(f, (x, -1, 1)).subs({x3:-x1,x4:-x2,x5:0,w3:w1,w4:w2}) I2 = (w1*(f.subs(x, x1)) + w2*(f.subs(x, x2)) + w3*(f.subs(x, x3)) + w4*(f.subs(x, x4)) + w5*(f.subs(x, x5))).subs({x3:-x1,x4:-x2,x5:0,w3:w1,w4:w2}) display("I1: ",I1) display("I2: ",I2) Eq1 = sp.expand(I1 - I2).coeff(a0) Eq2 = sp.expand(I1 - I2).coeff(a1) Eq3 = sp.expand(I1 - I2).coeff(a2) Eq4 = sp.expand(I1 - I2).coeff(a3) Eq5 = sp.expand(I1 - I2).coeff(a4) Eq6 = sp.expand(I1 - I2).coeff(a5) Eq7 = sp.expand(I1 - I2).coeff(a6) Eq8 = sp.expand(I1 - I2).coeff(a7) Eq9 = sp.expand(I1 - I2).coeff(a8) Eq10 = sp.expand(I1 - I2).coeff(a9) display("Eq1: ",Eq1) display("Eq2: ",Eq2) display("Eq3: ",Eq3) display("Eq4: ",Eq4) display("Eq5: ",Eq5) display("Eq6: ",Eq6) display("Eq7: ",Eq7) display("Eq8: ",Eq8) display("Eq9: ",Eq9) display("Eq10: ",Eq10) sol = list(sp.nonlinsolve([Eq1, Eq3, Eq5, Eq7, Eq9], [x1, x2, w1, w2, w5])) display(sol)

#### Eleventh-Degree Polynomial: Six Integration Points

For an eleventh-degree polynomial, , therefore, , so, 6 integration points are sufficient to exactly integrate an eleventh-degree polynomial. Following the procedure outlined above, the integration points along with the associated weight:

So, for functions that are very close to being eleventh-order polynomials, the following numerical integration scheme with 6 integration points can be employed:

The following is the Mathematica code used to find the above integration points and corresponding weights. Note that symmetry was employed to reduce the number of equations.

View Mathematica Codef = a0 + a1*x + a2*x^2 + a3*x^3 + a4*x^4 + a5*x^5 + a6*x^6 + a7*x^7 + a8*x^8 + a9*x^9 + a10*x^10 + a11*x^11; I1 = Integrate[f, {x, -1, 1}] I2 = w1*(f /. x -> x1) + w2*(f /. x -> x2) + w3*(f /. x -> x3) + w3*(f /. x -> -x3) + w2*(f /. x -> -x2) + w1*(f /. x -> -x1) Eq1 = Coefficient[I1 - I2, a0] Eq3 = Coefficient[I1 - I2, a2] Eq5 = Coefficient[I1 - I2, a4] Eq7 = Coefficient[I1 - I2, a6] Eq9 = Coefficient[I1 - I2, a8] Eq11 = Coefficient[I1 - I2, a10] a = Solve[{Eq1 == 0, Eq3 == 0, Eq5 == 0, Eq7 == 0, Eq9 == 0, Eq11 == 0}, {x1, x2, x3, w1, w2, w3}] Chop[N[a[[1]]]]

import numpy as np import sympy as sp from scipy.optimize import fsolve sp.init_printing(use_latex=True) x, x1, x2, x3 = sp.symbols('x x1 x2 x3') a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11 = sp.symbols('a0 a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 a11') w1, w2, w3, w4, w5 = sp.symbols('w1 w2 w3 w4 w5') f = a0 + a1*x + a2*x**2 + a3*x**3 + a4*x**4 + a5*x**5 + a6*x**6 + a7*x**7 + a8*x**8 + a9*x**9 + a10*x**10 + a11*x**11 I1 = sp.integrate(f, (x, -1, 1)) I2 = w1*(f.subs(x, x1)) + w2*(f.subs(x, x2)) + w3*(f.subs(x, x3)) + w3*(f.subs(x, -x3)) + w2*(f.subs(x, -x2)) + w1*(f.subs(x, -x1)) display("I1: ",I1) display("I2: ",I2) Eq1 = sp.expand(I1 - I2).coeff(a0) Eq2 = sp.expand(I1 - I2).coeff(a1) Eq3 = sp.expand(I1 - I2).coeff(a2) Eq4 = sp.expand(I1 - I2).coeff(a3) Eq5 = sp.expand(I1 - I2).coeff(a4) Eq6 = sp.expand(I1 - I2).coeff(a5) Eq7 = sp.expand(I1 - I2).coeff(a6) Eq8 = sp.expand(I1 - I2).coeff(a7) Eq9 = sp.expand(I1 - I2).coeff(a8) Eq10 = sp.expand(I1 - I2).coeff(a9) Eq11 = sp.expand(I1 - I2).coeff(a10) display("Eq1: ",Eq1) display("Eq2: ",Eq2) display("Eq3: ",Eq3) display("Eq4: ",Eq4) display("Eq5: ",Eq5) display("Eq6: ",Eq6) display("Eq7: ",Eq7) display("Eq8: ",Eq8) display("Eq9: ",Eq9) display("Eq10: ",Eq10) display("Eq11: ",Eq11) # Convert Sympy equations to a Python function using lambdify Equations = sp.lambdify(([x1, x2, x3, w1, w2, w3],), [Eq1, Eq3, Eq5, Eq7, Eq9, Eq11]) # Solve system of equations using Scipy's fsolve sol = fsolve(Equations, [-0.5,-1,-0.3,0,0,0]) display(sol)

We can continue determining more integration points (abscissae) and weights for Gauss quadrature. In this page, the integration points and weights up to 64 points are available. Note that 64 integration points of Gauss quadrature is associated with a polynomial of degree 127.

#### Implementation of Gauss Quadrature

Gauss quadrature is very easy to implement and provides very accurate results with very few computations. However, one drawback is that it is not applicable to data obtained experimentally as the values of the function at the specific integration points would not be necessarily available. In order to implement the Gauss integration scheme in Mathematica, first, the following table is created which contains the integration points and the corresponding weights. Row number contains the Gauss integration scheme with integration points. Then, a procedure is created in Mathematica whose input is a function , and the requested number of integration points. The procedure then calls the appropriate row in the “GaussTable” and calculates the weighted sum of the function evaluated at the appropriate integration points.

Clear[f, x] GaussTable = {{{"Integration Points"}, {"Corresponding Weights"}}, {{0}, {2}}, {{-1/Sqrt[3], 1/Sqrt[3]}, {1, 1}}, {{-Sqrt[3/5], 0, Sqrt[3/5]}, {5/9, 8/9, 5/9}}, {{-0.861136, -0.339981, 0.339981, 0.861136}, {0.347855, 0.652145, 0.652145, 0.347855}}, {{-0.90618, -0.538469, 0, 0.538469, 0.90618}, {0.236927, 0.478629, 0.568889, 0.478629, 0.236927}}, {{-0.93247, -0.661209, -0.238619, 0.238619, 0.661209, 0.93247}, {0.171324, 0.360762, 0.467914, 0.467914, 0.360762, 0.171324}}}; GaussTable // MatrixForm IG[f_, n_] := Sum[GaussTable[[n + 1, 2, i]]*f[GaussTable[[n + 1, 1, i]]], {i, 1, n}] f[x_] := x^9 + x^8 IG[f, 5.0]

import numpy as np import pandas as pd GaussTable = [[[0], [2]], [[-1/np.sqrt(3), 1/np.sqrt(3)], [1, 1]], [[-np.sqrt(3/5), 0, np.sqrt(3/5)], [5/9, 8/9, 5/9]], [[-0.861136, -0.339981, 0.339981, 0.861136], [0.347855, 0.652145, 0.652145, 0.347855]], [[-0.90618, -0.538469, 0, 0.538469, 0.90618], [0.236927, 0.478629, 0.568889, 0.478629, 0.236927]], [[-0.93247, -0.661209, -0.238619, 0.238619, 0.661209, 0.93247], [0.171324, 0.360762, 0.467914, 0.467914, 0.360762, 0.171324]]] display(pd.DataFrame(GaussTable, columns=["Integration Points", "Corresponding Weights"])) def IG(f, n): n = int(n) return sum([GaussTable[n - 1][1][i]*f(GaussTable[n - 1][0][i]) for i in range(n)]) def f(x): return x**9 + x**8 IG(f, 5.0)

##### Example

Calculate the exact integral of on the interval and find the relative error if a Gauss 1, 2, 3, and 4 integration points scheme is used. Also, calculate the number of rectangles required by the midpoint rule of the rectangle method to produce an error similar to that produced by the 3 point integration scheme, then evaluate the integral using the midpoint rule of the rectangle method with the obtained number of rectangles.

##### Solution

The exact integral is given by:

Employing a Gauss integration scheme with one integration point yields:

Employing a Gauss integration scheme with two integration points yields:

Employing a Gauss integration scheme with three integration points yields:

Employing a Gauss integration scheme with four integration points yields:

Using a three point integration scheme yielded an error of:

Equating the upper bound for the error in the midpoint rule to the error produced by a three-point Gauss integration scheme yields:

With , and . Therefore:

Therefore, the required number of rectangles is:

With 75 rectangles the rectangle method yields an integral of . It is important to note how accurate the Gauss integration scheme is compared to the rectangle method. For this example, it would take 75 divisions (equivalent to 75 integration points) to reach an accuracy similar to the Gauss integration scheme with only 3 points!

View Mathematica CodeClear[f, x, I2] GaussTable = {{{"Integration Points"}, {"Corresponding Weights"}}, {{0}, {2}}, {{-1/Sqrt[3], 1/Sqrt[3]}, {1, 1}}, {{-Sqrt[3/5], 0, Sqrt[3/5]}, {5/9, 8/9, 5/9}}, {{-0.861136, -0.339981, 0.339981, 0.861136}, {0.347855, 0.652145, 0.652145, 0.347855}}, {{-0.90618, -0.538469, 0, 0.538469, 0.90618}, {0.236927, 0.478629, 0.568889, 0.478629, 0.236927}}, {{-0.93247, -0.661209, -0.238619, 0.238619, 0.661209, 0.93247}, {0.171324, 0.360762, 0.467914, 0.467914, 0.360762, 0.171324}}}; GaussTable // MatrixForm IG[f_, n_] := Sum[GaussTable[[n + 1, 2, i]]*f[GaussTable[[n + 1, 1, i]]], {i, 1, n}]; f[x_] := Cos[x] Iexact = Integrate[Cos[x], {x, -1, 1.0}] Itable = Table[{i, N[IG[f, i]], (Iexact - IG[f, i])/Iexact}, {i, 1, 6}]; Title = {"Number of Integration Points", "Numerical Integration Results", "Relative Error"}; Itable = Prepend[Itable, Title]; Itable // MatrixForm h = Sqrt[0.00006*24/2/1] n = Ceiling[2/h] I2[f_, a_, b_, n_] := (h = (b - a)/n; Sum[f[a + (i - 1/2)*h]*h, {i, 1, n}]) I2[f, -1.0, 1, n]

import numpy as np import sympy as sp import pandas as pd from scipy import integrate GaussTable = [[[0], [2]], [[-1/np.sqrt(3), 1/np.sqrt(3)], [1, 1]], [[-np.sqrt(3/5), 0, np.sqrt(3/5)], [5/9, 8/9, 5/9]], [[-0.861136, -0.339981, 0.339981, 0.861136], [0.347855, 0.652145, 0.652145, 0.347855]], [[-0.90618, -0.538469, 0, 0.538469, 0.90618], [0.236927, 0.478629, 0.568889, 0.478629, 0.236927]], [[-0.93247, -0.661209, -0.238619, 0.238619, 0.661209, 0.93247], [0.171324, 0.360762, 0.467914, 0.467914, 0.360762, 0.171324]]] display(pd.DataFrame(GaussTable, columns=["Integration Points", "Corresponding Weights"])) def IG(f, n): n = int(n) return sum([GaussTable[n - 1][1][i]*f(GaussTable[n - 1][0][i]) for i in range(n)]) def f(x): return np.cos(x) Iexact, error = integrate.quad(f, -1, 1) print("Iexact: ",Iexact) Itable = [[i + 1, sp.N(IG(f, i + 1)), (Iexact - IG(f, i + 1))/Iexact] for i in range(6)] Itable = pd.DataFrame(Itable, columns=["Number of Integration Points", "Numerical Integration Results", "Relative Error"]) display(Itable) h = np.sqrt(0.00006*24/2/1) print("h: ",h) n = np.ceil(2/h) print("n: ",n) def I2(f, a, b, n): h = (b - a)/n return sum([f(a + (i + 1/2)*h)*h for i in range(int(n))]) I2(f, -1.0, 1, n)

##### Alternate Integration Limits

The Gauss integration scheme presented above applies to functions that are integrated over the interval . A simple change of variables can be used to allow integrating a function where . In this case, the linear relationship between and can be expressed as:

Therefore:

Which means the integration can be converted from integrating over to integrating over as follows:

where

Therefore:

The integration can then proceed as per the weights and values of the integration points with given by:

##### Example

Calculate the exact integral of on the interval and find the absolute relative error if a Gauss 1, 2, 3, and 4 integration points scheme is used.

##### Solution

First, to differentiate between the given function limits, and the limits after changing variables, we will assume that the function is given in terms of as follows:

The exact integral is given by:

Using a linear change of variables we have:

Therefore:

Employing a Gauss integration scheme with one integration point yields:

Employing a Gauss integration scheme with two integration points yields:

Employing a Gauss integration scheme with three integration points yields:

where,

In the following Mathematica code, the procedure IGAL applies the Gauss integration scheme with integration points to a function on the interval

View Mathematica CodeGaussTable = {{{"Integration Points"}, {"Corresponding Weights"}}, {{0}, {2}}, {{-1/Sqrt[3], 1/Sqrt[3]}, {1, 1}}, {{-Sqrt[3/5], 0, Sqrt[3/5]}, {5/9, 8/9, 5/9}}, {{-0.861136, -0.339981, 0.339981, 0.861136}, {0.347855, 0.652145, 0.652145, 0.347855}}, {{-0.90618, -0.538469, 0, 0.538469, 0.90618}, {0.236927, 0.478629, 0.568889, 0.478629, 0.236927}}, {{-0.93247, -0.661209, -0.238619, 0.238619, 0.661209, 0.93247}, {0.171324, 0.360762, 0.467914, 0.467914, 0.360762, 0.171324}}}; GaussTable // MatrixForm IGAL[f_, n_, a_, b_] := Sum[(b - a)/2*GaussTable[[n + 1, 2, i]]*f[(b - a)/2*(GaussTable[[n + 1, 1, i]] + 1) + a], {i, 1, n}] f[x_] := x*E^x Iexact = Integrate[f[x], {x, 0., 3}] Itable = Table[{i, N[IGAL[f, i, 0, 3]], (Iexact - IGAL[f, i, 0, 3])/Iexact}, {i, 1, 6}]; Title = {"Number of Integration Points", "Numerical Integration Results", "Relative Error"}; Itable = Prepend[Itable, Title]; Itable // MatrixForm

import numpy as np import sympy as sp import pandas as pd from scipy import integrate GaussTable = [[[0], [2]], [[-1/np.sqrt(3), 1/np.sqrt(3)], [1, 1]], [[-np.sqrt(3/5), 0, np.sqrt(3/5)], [5/9, 8/9, 5/9]], [[-0.861136, -0.339981, 0.339981, 0.861136], [0.347855, 0.652145, 0.652145, 0.347855]], [[-0.90618, -0.538469, 0, 0.538469, 0.90618], [0.236927, 0.478629, 0.568889, 0.478629, 0.236927]], [[-0.93247, -0.661209, -0.238619, 0.238619, 0.661209, 0.93247], [0.171324, 0.360762, 0.467914, 0.467914, 0.360762, 0.171324]]] display(pd.DataFrame(GaussTable, columns=["Integration Points", "Corresponding Weights"])) def IGAL(f, n, a, b): n = int(n) return sum([(b - a)/2*GaussTable[n - 1][1][i]*f((b - a)/2*(GaussTable[n - 1][0][i] + 1) + a) for i in range(n)]) def f(x): return x*np.exp(x) Iexact, error = integrate.quad(f, 0, 3) print("Iexact: ",Iexact) Itable = [[i + 1, sp.N(IGAL(f, i + 1, 0, 3)), (Iexact - IGAL(f, i + 1, 0, 3))/Iexact] for i in range(6)] Itable = pd.DataFrame(Itable, columns=["Number of Integration Points", "Numerical Integration Results", "Relative Error"]) display(Itable)

Hello! Very helpful article! I do have one question though. For the python code, what if I would like to solve a double integral? As an example, what if I have r^2*cos(theta)drdtheta and I would like to integrate in terms of both r and theta (from 0 to 1 and 0 to pi respectively)? I do see how one can do this in the Mathematica version (likely due to my greater experience using Mathematica) but I keep running into issues attempting this in Python. Any help would be greatly appreciated

You can post your code here (or email it to me) and I can take a look.

It should be straightforward in any language. You can review the following page:

https://engcourses-uofa.ca/books/introduction-to-solid-mechanics/finite-element-analysis/one-and-two-dimensional-isoparametric-elements-and-gauss-integration/gauss-integration/#gauss-integration-over-two-dimensional-domains-6

Here is the code I came up with. I’m assuming one of the problems is how I define “theta” (if you run “ampf” it’s clear that the function is not actually integrating in terms of theta). I have slowly learned that defining variables in python is a bit more involved than mathematica!

import numpy as np

from matplotlib import pyplot as plt

import matplotlib as mpl

import sympy as sp

from mpmath import *

import pandas as pd

from scipy import integrate

GaussTable = [[[0], [2]], [[-1/np.sqrt(3), 1/np.sqrt(3)], [1, 1]], [[-np.sqrt(3/5), 0, np.sqrt(3/5)], [5/9, 8/9, 5/9]], [[-0.861136, -0.339981, 0.339981, 0.861136], [0.347855, 0.652145, 0.652145, 0.347855]], [[-0.90618, -0.538469, 0, 0.538469, 0.90618], [0.236927, 0.478629, 0.568889, 0.478629, 0.236927]], [[-0.93247, -0.661209, -0.238619, 0.238619, 0.661209, 0.93247], [0.171324, 0.360762, 0.467914, 0.467914, 0.360762, 0.171324]]]

display(pd.DataFrame(GaussTable, columns=[“Integration Points”, “Corresponding Weights”]))

def IG(f, n):

n = int(n)

return sum([GaussTable[n – 1][1][i]*f(GaussTable[n – 1][0][i]) for i in range(n)])

theta = sp.symbols(“theta”)

def f(r): return r**2*sp.cos(theta)

n = 200

def I2(f, a, b, n):

h = (b – a)/n

return sum([f(a + (i + 1/2)*h)*h for i in range(int(n))])

ampr=I2(f, 0, 1, n)

def g(theta): return ampr

def I2(g, a, b, n):

h = (b – a)/n

return sum([g(a + (i + 1/2)*h)*h for i in range(int(n))])

ampf=I2(g,0,3.14159265358,n)

The integration function should look different. Please check this app that I created for this:

https://mecsimcalc.com/app/0505681/gauss_integration

You can find the code under: “app doc”. It seems to work properly up to th2=1.0Pi. I tried th2=2.0Pi and it didn’t give the expected output and I might have made a mistake. If you find a mistake, let me know.

Thanks for the code! I played around with it and got something that seems to mostly work. The only issue I have found is that I get an error when n1 and n2 are greater than 6. Not fully sure why this is the case. For the example I give below n1=n2=6 is fine but I imagine I would need more steps for more complicated integrals. Any thoughts?

import numpy as np

from matplotlib import pyplot as plt

import matplotlib as mpl

import sympy as sp

from mpmath import *

import pandas as pd

from scipy import integrate

GaussTable = [[[0], [2]], [[-1/np.sqrt(3), 1/np.sqrt(3)], [1, 1]], [[-np.sqrt(3/5), 0, np.sqrt(3/5)], [5/9, 8/9, 5/9]], [[-0.861136, -0.339981, 0.339981, 0.861136], [0.347855, 0.652145, 0.652145, 0.347855]], [[-0.90618, -0.538469, 0, 0.538469, 0.90618], [0.236927, 0.478629, 0.568889, 0.478629, 0.236927]], [[-0.93247, -0.661209, -0.238619, 0.238619, 0.661209, 0.93247], [0.171324, 0.360762, 0.467914, 0.467914, 0.360762, 0.171324]]]

display(pd.DataFrame(GaussTable, columns=[“Integration Points”, “Corresponding Weights”]))

def IGAL(f, n1, r1, r2, n2,th1,th2):

n1 = int(n1)

n2 = int(n2)

return sum([sum([(r2 – r1)/2*GaussTable[n1 – 1][1][i]*(th2 – th1)/2*GaussTable[n2 – 1][1][j]*f((r2 – r1)/2*(GaussTable[n1 – 1][0][i] + 1) + r1,(th2 – th1)/2*(GaussTable[n2 – 1][0][j] + 1) + th1) for i in range(n1)]) for j in range(n2)])

def f(r,th):

return r**2*np.sin(th)

n1 = 6

n2 = 6

ampf=IGAL(f,n1,0,1,n2,0,np.pi)

print(ampf)

The Gauss Table data provided in the code is only up to n=6

Ohhh I see! Looks like I’m going to have to figure out how to calculate a much larger Gauss Table going forward. For the integration I plan to do, I may likely need n=200! (*not a factorial I’m just excited haha!)