11. Numerical Integration

11.1 Overview

Let \(f(x)\) be a function defined over the interval \([a, b]\), and consider its definite integral:
\[
I = \int_a^b f(x) dx
\] to be computed numerically.

There are two major types of numerical integration programs. One is the function-input type, which is used when the value of the integrand can be calculated at any point within the interval upon request. The other is the data-input type, which is used when the values of the integrand are known only at a few points within the interval.

11.2 Methods for Numerical Integration

Example Problem

Compute the following definite integral:
\[
I = \int_0^4 1/(x + 1) dx
\] The function shape of \(1/(x + 1)\) is as follows:

This function can be analytically integrated, yielding \(ln(x + 1) + C\) (where \(C\) is the integration constant). Evaluating for \(a = 0\) and \(b = 4\) gives \(I = 1.609437912\).

11.2.1 Function Input Type

There are many formulas for approximating integrals, such as trapezoidal rule and Simpson’s rule.

11.2.1.1 Trapezoidal Rule

Approximates the integral using a linear interpolation between \(x = a\) and \(x = b\), treating the trapezoid area as an approximation:
\[
I = (f(a) + f(b))(b – a)/2
\] Applying it to the example gives:

The computed result is \(I = 2.4000\) (error: 49.1%).

11.2.1.2 Simpson’s Rule

Evaluates the function at an additional midpoint, performing quadratic interpolation using three points:
\[
I = (f(a) + 4f((a + b)/2) + f(b))(b – a)/6
\] Applying it to the example:

The computed result is \(I = 1.6889\) (error: 4.94%), significantly more accurate than the trapezoidal rule. However, Simpson’s Rule requires three function evaluations, while the trapezoidal rule only requires two.

The computational cost of numerical integration (1/(computation speed)) can generally be considered to depend on the number of evaluations of \(f(x)\). In this case, Simpson’s rule requires 1.5 times the computational cost compared to the trapezoidal rule, but it provides more than 10 times the accuracy.

11.2.1.3 Composite Rules

To improve accuracy, the integral interval can be subdivided into multiple smaller intervals, applying the trapezoidal or Simpson’s rule to each subinterval. The length of each subinterval may be different. However, in this example, the equally spaced subintervals with the length \(h\) are assumed.

Let the interval \([a, b]\) be divided into \(n\) subintervals \([x_i, x_{i+1}] (i = 1 \sim n)\):
\[
I = \sum_{i=1}^n I_i
\] where,
\[
I_i = \int_{x_i}^{x_{i+1}}f(x) dx \\
x_1 = a, x_{n+1} = b, h = (b – a)/n
\] Applying the trapezoidal rule to each subinterval:
\[
I_i = h(f(x_i + f(x_{i+1}))/2
\] This is called the composite trapezoidal rule. However, the term “trapezoidal rule” often refers to the composite trapezoidal rule.

Dividing the example interval into four subintervals \((m = 4)\) and applying the trapezoidal rule:

The computed result is \(I = 1.6833\) (error 4.59%), which is about 10 times better than the previous error of 49.1% before partitioning. The number of function evaluations for \(f(x)\) increased to 5, which is 2.5 times more.

Similarly, applying Simpson’s rule to subintervals:
\[
I_i = h(f(x_i + 4f((x_i + x_{i+1})/2) + f(x_{i+1}))/6
\] This is called the composite Simpson’s rule. However, the term “Simpson’s rule” often refers to the composite Simpson’s rule.

Applying the example with a two-part division \((m = 2)\) so that the number of function evaluations for \(f(x)\) matches the composite trapezoidal rule (5 evaluations), the result is as follows.

The computed result is \(I = 1.6222\) (error 0.79%), which is approximately five times better than the previous error of 4.94% before partitioning. Even with the same number of function evaluations (half the number of partitions), Simpson’s rule demonstrates significantly higher accuracy compared to the trapezoidal rule.

11.2.1.4 Gaussian Quadrature

Gaussian quadrature (also known as Gauss-Legendre quadrature) offers higher accuracy than Simpson’s rule.

This method interpolates the integral using Legendre polynomials, which are one of the orthogonal polynomials.

In this formula, evaluation points \(x_i\) are fixed at the zeros of the Legendre polynomial, and the integral is approximated using weights \(w_i\):
\[
I = ((b – a)/2)\sum_{i=1}^n w_if((a + b)/2 + ((b – a)/2)x_i)
\] Values of \(x_i\) and \(w_i\) can be found in textbooks or formula tables. The following are values for n = 2 and n = 3, referred to as 2-point Gaussian quadrature and 3-point Gaussian quadrature:
\[
\begin{array}{ccccccc}
n & x1 & x2 & x3 & w1 & w2 & w3 \\ \hline
2 & -\sqrt{1/3} & \sqrt{1/3)} & – & 1 & 1 & – \\
3 & -\sqrt{3/5} & 0 & \sqrt{3/5} & 5/9 & 8/9 & 5/9 \\
\hline
\end{array}
\] Note that \(x_i\) does not include the interval endpoints \(a\) and \(b\).

Applying 2-point Gaussian quadrature to the example:

The number of function evaluations is 2 with the error 2.75%.

Applying 3-point Gaussian quadrature:

The number of function evaluations is 3 with the error 0.42%.

Gaussian quadrature exhibits significantly higher precision than the trapezoidal and Simpson’s rules. Further tests for \(n = 4\) and \(5\) show errors 0.063% and 0.009%, respectively.

In the case of Gaussian quadrature, to improve accuracy, either \(n\) should be increased, or the integration interval should be divided into subintervals and the quadrature formula should be applied to each subinterval.

11.2.1.5 Adaptive Quadrature Methods

If integration error estimates are available, computations can be optimized by dynamically adjusting the number of subintervals, ensuring the required accuracy while minimizing computation cost. This approach is known as adaptive quadrature.

To estimate the error, one possible approach is to vary the number of subintervals and compare the results. When using the trapezoidal rule or Simpson’s rule with equally spaced subintervals, doubling the number of subintervals (halving the width of each) allows for the reuse of previously computed \(f(x)\) values at half of the points, thereby reducing computational cost.

In contrast, with the Gaussian quadrature formula, changing \(n\) requires recalculating \(f(x)\) at all points, making it computationally inefficient. To address this, the Gauss-Kronrod quadrature formula was proposed, allowing two different \(n\) values to be obtained within a single formula. The approach is as follows.
\[
I = ((b – a)/2)(\sum_{i=1}^n a_if((a + b)/2 + ((b – a)/2)x_i) + \sum_{j=1}^{n+1} b_jf((a + b)/2 + ((b – a)/2)y_j))
\] where \(a_i\) and \(b_j\) are Gauss-Kronrod quadrature coefficients. Do not confuse with the interval endpoints \(a\) and \(b\).

The \(x_i\) values are the same as those in the Gaussian quadrature formula. Therefore, by combining them with the weights \(w_i\) from the Gaussian quadrature, the integral value using the \(n\)-point Gaussian quadrature (denoted as \(G_n\)) can be obtained. Note that the values of \(a_i\) and \(w_i\) are different.

The integral value using the \((2n + 1)\)-point Gauss-Kronrod quadrature (denoted as \(K_{2n+1}\)) can be computed by the full formula.

\(K_{2n+1}\) has lower accuracy compared to \(G_{2n+1}\) despite having the same number of function evaluations for \(f(x)\). In other words, the Gauss-Kronrod quadrature formula introduces an error estimation feature at the cost of slightly reduced accuracy compared to the Gaussian quadrature formula.

For \(n = 2\) (5-point Gauss-Kronrod quadrature), coefficients are:
\[
\begin{array}{cccc}
\hline
x_i & -0.5773502691896258 & 0.5773502691896258 \\
a_i & 0.4909090909090909 & 0.4909090909090909 \\
y_j & -0.9258200997725515 & 0 & 0.9258200997725515 \\
b_j & 0.1979797979797980 & 0.6222222222222222 & 0.1979797979797980 \\
\hline
\end{array}
\] From \(G_n\) and \(K_{2n+1}\), the integration error can be estimated.

Using 5-point Gauss-Kronrod quadrature for the example resulted in 0.012% error. It is an intermediate precision between 4-point and 5-point Gaussian quadrature.

Nowadays, numerical integration commonly employs adaptive quadrature routines to ensure a certain level of accuracy. In particular, programs that apply the Gauss-Kronrod quadrature formula to each subinterval for error estimation and optimize the width of each subinterval to reduce computational cost while maintaining accuracy are widely used.

11.2.2 Data Input Type

The integral value can be obtained by interpolating the given data points and calculating the integral of the resulting interpolation function.

There are methods that use polynomial interpolation and those that use spline interpolation. Accuracy varies depending on input data quality, but this method typically does not yield very high precision results.

11.3 Numerical Integration Using XLPack

The VBA subroutine Qk15 performs numerical integration using the 15-point Gauss-Kronrod quadrature formula, calculating the integral value and estimated error in 15 function evaluations. The VBA subroutine Qag is an adaptive quadrature routine based on the Gauss-Kronrod formula. If the function is smooth, Qk15 provides sufficient accuracy. However, if the function is not smooth or higher precision is required with more function evaluations, Qag is recommended.

For infinite-range integration, the VBA subroutine Qagi can be used. It also applies the adaptive Gauss-Kronrod quadrature formula, transforming semi-infinite integrals into finite interval integrals over [0, 1].

Both Qag and Qagi are accessible via the XLPack solver.

By combining the VBA subroutine Pchse or the worksheet function WPchse (see “5. Interpolation Methods”) with the VBA subroutine Pchia or the worksheet function WPchia, it can be used as a numerical integration program for data input type calculations.

Below, examples are provided for solving the previous integral problem.

11.3.1 Solving Using VBA Program (Function Input) (1)

An example using Qag and Qk15:

Function F(X As Double) As Double
    F = 1 / (1 + X)
End Function

Sub Start()
    Dim A As Double, B As Double, S As Double, Info As Long
    Dim AbsErr As Double, Neval As Long
    '--- Input data
    A = Cells(4, 2): B = Cells(5, 2)
    '--- Compute integration by Qag
    Call Qag(AddressOf F, A, B, S, Info, AbsErr, Neval)
    Cells(11, 2) = Info
    If Info = 0 Then
        Cells(8, 2) = S
        Cells(9, 2) = AbsErr
        Cells(10, 2) = Neval
    End If
    '--- Compute integration by Qk15
    Call Qk15(AddressOf F, A, B, S, AbsErr)
    Cells(8, 3) = S
    Cells(9, 3) = AbsErr
End Sub

An external function defines the integrand, while Qag and Qk15 compute the integral using the specified interval.

Executing this program gives the following results:

Qag produced a result with an estimated error of \(7.38 \times 10^{-13}\) after 75 function calls (internally calling Qk15 five times). Qk15 had an estimated error of \(2.07 \times 10^{-5}\) after 15 function calls. The estimated errors were calculated conservatively, and the actual errors were \(2.22 \times 10^{-16}\) and \(1.20 \times 10^{-11}\), respectively.

11.3.2 Solving Using VBA Program (Function Input) (2)

An example using Qag_r and Qk15_r (reverse communication interface):

Sub Start()
    Dim A As Double, B As Double, S As Double, Info As Long
    Dim AbsErr As Double, Neval As Long
    Dim XX As Double, YY As Double, IRev As Long
    '--- Input data
    A = Cells(4, 2): B = Cells(5, 2)
    '--- Compute integration by Qag_r
    IRev = 0
    Do
        Call Qag_r(A, B, S, Info, XX, YY, IRev, AbsErr, Neval)
        If IRev = 1 Then YY = 1 / (1 + XX ^ 2)
    Loop While IRev <> 0
    Cells(11, 2) = Info
    If Info = 0 Then
        Cells(8, 2) = S
        Cells(9, 2) = AbsErr
        Cells(10, 2) = Neval
    End If
    '--- Compute integration by Qk15_r
    IRev = 0
    Do
        Call Qk15_r(A, B, S, XX, YY, IRev, AbsErr)
        If IRev = 1 Then YY = 1 / (1 + XX ^ 2)
    Loop While IRev <> 0
    Cells(8, 3) = S
    Cells(9, 3) = AbsErr
End Sub

Here, the function values are computed inside the loop when IRev = 1, and then fed back into Qag_r or Qk15_r. More details on RCI are available here.

Executing this program produces the same results.

11.3.3 Solving Using Worksheet Functions (Data Input)

For data-based numerical integration, spline interpolation is used.

In the same example as above, the functional form is unknown, and only the function values at \(x = 0, 0.5, 1, \dots, 4\) are given.
\[
\begin{array}{cccc}
\hline
x & f(x) \\ \hline
0 & 1 \\
0.5 & 0.666667 \\
1 & 0.5 \\
1.5 & 0.4 \\
2 & 0.333333 \\
2.5 & 0.285714 \\
3 & 0.25 \\
3.5 & 0.222222 \\
4 & 0.2 \\
\hline
\end{array}
\] First, use the WPchse function to determine the spline function coefficients (see “5. Interpolation”).

Once the coefficients have been obtained, the WPchia function can be used to calculate the integral value. The required parameters for WPchia are A, B, N, X, Y, and D. A and B represent the integration range (0 to 4 in this example). N is the number of data points (9 in this example). X and Y define the cell range for the interpolation data, while D specifies the cell range for the spline coefficients.

Executing this function yields 1.61082…, achieving three-digit accuracy from nine data points.

Note – The example worksheet also includes an integration example using spline interpolation with VBA.

11.3.4 Solving Using the Solver

You can also solve it using the “Quadrature (Finite Interval)” feature of the XLPack Solver Add-in. The integrand function (=1/(1+B3^2)) is entered in cell B4.

For more information about the solver, refer to this link.