Abstract

This work presents a matrix-based reformulation of the composite Simpson’s 1/3 rule for the numerical evaluation of double integrals over rectangular domains with continuous integrands and prescribed non-uniformly spaced nodes. The classical formulation is usually expressed through nested summations or through a matrix of order $\left(n+1\right)\times \left(m+1\right)$ , whose size increases with the number of discretization nodes. The proposed approach reorganizes the contributions of the subrectangles into a fixed 3 × 3 matrix whose nine entries are defined by weighted double summations over the computational mesh. Consequently, the numerical approximation can be written compactly as the product of matrices of orders 1 × 3, 3 × 3, and 3 × 1, with the vectors containing the Simpson coefficients $\left(1,4,1\right)$ . The fixed-size representation provides a compact organization of the integration procedure and can reduce the auxiliary memory required for the matrix representation, while the asymptotic computational complexity remains $O\left(nm\right)$ . Numerical results for non-uniform meshes show that the formulation provides accurate approximations consistent with the behavior of the composite Simpson’s 1/3 rule. The proposed method is intended as a simple and flexible alternative for double numerical integration when data are available on prescribed non-uniform rectangular meshes.

Keywords

Matrix Formulation of Numerical Integration, Double Numerical Integration, Composite Simpson’s 1/3 Rule

Share and Cite:

1. Introduction

Otto and Denier [1] state that, although mathematics is an elegant and precise discipline, there are many problems that simply do not have analytical solutions, others whose exact solutions are unknown, and still others whose analytical solutions are too cumbersome to obtain. In such situations, numerical analysis provides practical techniques that transform these problems into finite-dimensional ones that can be solved computationally. Likewise, Süli and Mayers [2] point out that “The problem of evaluating definite integrals arises both in mathematics and beyond it, in many areas of science and engineering”. However, many definite integrals cannot be evaluated analytically, making numerical integration an essential tool in scientific computing. With respect to the best-known numerical integration methods, Nieves and Domínguez [3] state that “Commonly used integration methods can be classified into two groups: those that use given values of the function f(x) at equally spaced abscissas, known as Newton-Cotes formulas, and those that use values of f(x) at unequally spaced abscissas determined by specific properties of families of orthogonal polynomials, known as Gaussian quadrature formulas”. Although Grant and Talvila [4] observe that “There is a vast literature on numerical integration of functions of a single variable; however, the literature on the numerical integration of functions of several variables is sparse”, considerable research has been devoted to the development of efficient numerical integration techniques for multidimensional problems. Among the earliest contributions, Lee et al. [5] proposed an optimization-based approach for the numerical evaluation of two-dimensional integrals by formulating the construction of quadrature rules as an optimization problem. More recently, Keshavarzzadeh et al. [6] developed a computational framework for constructing positive quadrature rules in multiple dimensions based on moment-matching and optimization techniques, generating efficient integration schemes for general geometries while preserving polynomial exactness. Another active line of research is based on wavelet approximations. Aziz et al. [5] developed quadrature rules for numerical integration using Haar wavelets and hybrid functions, while Dehda [7] proposed a method for the numerical approximation of double integrals based on Coifman wavelets, achieving smaller absolute errors than methods based on Haar wavelets, CAS wavelets, and hybrid functions. Although these methods provide accurate approximations, they rely on wavelet theory and operational matrices, resulting in formulations that are mathematically more involved. Recent studies have also focused on improving classical quadrature and cubature rules. In [8], a hybrid numerical integration method combining Simpson’s rule and Gaussian quadrature was proposed for single-variable integrals, demonstrating that the enhancement of classical quadrature schemes remains an active research topic. Extending this idea to double integrals, Mohanty and Alahmadi [9] proposed a generalized cubature rule by combining the Gauss-Legendre two-point cubature rule with Simpson’s 1/3 cubature rule, obtaining a higher-degree precision method while maintaining computational efficiency. In a different direction, Zhong and Feng [10] introduced an efficient tensor-product reformulation for quasi-Monte Carlo lattice rules to accelerate the numerical approximation of high-dimensional integrals through an algebraic reformulation. Unlike the aforementioned approaches, which develop optimization-based methods, wavelet-based formulations, hybrid quadrature or cubature rules, or tensor-product reformulations, the present work develops a simple matrix reformulation of the classical composite Simpson’s 1/3 rule for double numerical integration over rectangular domains with non-uniformly spaced data. The proposed formulation leads to a fixed 3 × 3 matrix representation whose size is independent of the number of discretization nodes, thereby providing a compact implementation of the classical composite Simpson’s rule.

The main contribution of this work is the development of a simple and practical method for the numerical evaluation of double integrals over rectangular domains with continuous integrands and non-uniformly spaced data. This method is based on novel matrix-based reformulation of the standard composite Simpson’s 1/3 rule, given by:

$\begin{array}{l}{\displaystyle {\int }_c^d{\displaystyle {\int }_a^{b}f\left(x,y\right)\text{d}x\text{d}y}}\\ \approx \frac{h_1{h}_2}{9}\left[f\left(x_0,y_0\right)+4{\displaystyle \sum _{\begin{array}{c}j=1\\ \text{Δ}j=2\end{array}}^{m-1}f}\left(x_0,y_j\right)+2{\displaystyle \sum _{\begin{array}{c}j=2\\ \text{Δ}j=2\end{array}}^{m-2}f}\left(x_0,y_j\right)+f\left(x_0,y_m\right)\right]\\ +\frac{4h_1{h}_2}{9}{\displaystyle \sum _{\begin{array}{c}i=1\\ \text{Δ}i=2\end{array}}^{n-1}\left[f\left(x_i,y_0\right)+4{\displaystyle \sum _{\begin{array}{c}j=1\\ \text{Δ}j=2\end{array}}^{m-1}f}\left(x_i,y_j\right)+2{\displaystyle \sum _{\begin{array}{c}j=2\\ \text{Δ}j=2\end{array}}^{m-2}f}\left(x_i,y_j\right)+f\left(x_i,y_m\right)\right]}\\ +\frac{2h_1{h}_2}{9}{\displaystyle \sum _{\begin{array}{c}i=2\\ \text{Δ}i=2\end{array}}^{n-2}\left[f\left(x_i,y_0\right)+4{\displaystyle \sum _{\begin{array}{c}j=1\\ \text{Δ}j=2\end{array}}^{m-1}f}\left(x_i,y_j\right)+2{\displaystyle \sum _{\begin{array}{c}j=2\\ \text{Δ}j=2\end{array}}^{m-2}f}\left(x_i,y_j\right)+f\left(x_i,y_m\right)\right]}\\ +\frac{h_1{h}_2}{9}\left[f\left(x_n,y_0\right)+4{\displaystyle \sum _{\begin{array}{c}j=1\\ \text{Δ}j=2\end{array}}^{m-1}f}\left(x_n,y_j\right)+2{\displaystyle \sum _{\begin{array}{c}j=2\\ \text{Δ}j=2\end{array}}^{m-2}f}\left(x_n,y_j\right)+f\left(x_n,y_m\right)\right]\end{array}$

A key novelty of the proposed approach lies in transforming the classical summation-based formulation into a compact matrix representation. Specifically, the double integral can be expressed as a structured matrix product, thereby providing a unified algebraic framework for its evaluation.

More precisely, the numerical approximation can be written in the form

${\displaystyle {\int }_c^d{\displaystyle {\int }_a^{b}f\left(x,y\right)\text{d}x\text{d}y}}\approx \left(\begin{array}{lll}1\hfill & 4\hfill & 1\hfill \end{array}\right)M\left(\begin{array}{c}1\\ 4\\ 1\end{array}\right)$

where

$M=\left(\begin{array}{ccc}{M}_{11}& M_{12}& M_{13}\\ M_{21}& M_{22}& M_{23}\\ M_{31}& M_{32}& M_{33}\end{array}\right)$

The coefficients 1, 4, 1 correspond to those of Simpson’s 1/3 numerical integration rule. The quantities $M_{jk}$ , for $j,k=1,2,3$ , represent double summations which, to the best of our knowledge, have not been previously reported in the literature. Consequently, the numerical evaluation of the definite double integral is equivalent to computing its associated matrix $M$ . Unlike the classical formulation, which involves nested summations over all nodes, the proposed method reduces the problem to operations involving a fixed 3 × 3 matrix, regardless of the number of discretization points. This feature not only simplifies implementation but also provides a more transparent interpretation in terms of linear algebra, making the method compact and conceptually insightful.

In section 2, we present the basic theory necessary as a starting point for the development of the main idea of this work.

Equation (11) of Section 2.4, a particular case of Equation (10) of Section 2.3, together with the concept of composite integration, constitutes the basis for the main contribution of this work presented in Section 3.

2. Theory

2.1. Simpson’s 1/3 Rule

Simpson’s 1/3 rule is based upon a quadratic interpolation of a single variable function and is a special case of closed Newton-Cotes formulas. It needs two equal subdivisions of the integration range $\left[a,b\right]$ , and hence the integrand is evaluated in three points called nodes $x_0,x_1$ , and $x_2$ , distributed symmetrically in the interval $\left[a,b\right]$ . Then, Simpson’s 1/3 rule is written as follows:

${\displaystyle {\int }_a^{b}f\left(x\right)\text{d}x}\approx \frac{h}{3}\left[f\left(x_0\right)+4f\left(x_1\right)+f\left(x_2\right)\right]$ (1)

where $h=\frac{b-a}{2}$ is the step size, $x_0=a$ and $x_2=b$ .

2.2. Matrix Formulation of the Composite Simpson’s 1/3 Integration Rule for Functions of One Variable

By using Equation (1) we can derive the formula for the composite integration Simpson 1/3, expressed in terms of sums and given by

${\displaystyle {\int }_a^{b}f\left(x\right)\text{d}x}\approx \frac{h}{3}\left[f\left(x_0\right)+4{\displaystyle \sum _{\begin{array}{c}i=1\\ \text{Δ}i=2\end{array}}^{n-1}f}\left(x_i\right)+2{\displaystyle \sum _{\begin{array}{c}i=2\\ \text{Δ}i=2\end{array}}^{n-2}f}\left(x_i\right)+f\left(x_n\right)\right]$ (2)

where $h=\frac{b-a}{n}$ is the step size.

Expanding the summations in Equation (2) into their explicit weighted form, for example when n = 6, we obtain

${\displaystyle {\int }_a^{b}f\left(x\right)\text{d}x}\approx \frac{h}{3}\left[f\left(x_0\right)+4f\left(x_1\right)+2f\left(x_2\right)+4f\left(x_3\right)+2f\left(x_4\right)+4f\left(x_5\right)+f\left(x_6\right)\right]$

This expression naturally leads to the following matrix formulation:

${\displaystyle {\int }_a^{b}f\left(x\right)\text{d}x}\approx \frac{h}{3}\left(\begin{array}{lllllll}1\hfill & 4\hfill & 2\hfill & 4\hfill & 2\hfill & 4\hfill & 1\hfill \end{array}\right)\left(\begin{array}{c}f\left(x_0\right)\\ f\left(x_1\right)\\ f\left(x_2\right)\\ f\left(x_3\right)\\ f\left(x_4\right)\\ f\left(x_5\right)\\ f\left(x_6\right)\end{array}\right)$ (3)

The coefficient vector has 1 as its first and last entries, while the intermediate entries alternate between 4 and 2, starting and ending with 4. By generalizing Equation (3), we obtain the following matrix formulation of Equation (2).

${\displaystyle {\int }_a^{b}f\left(x\right)\text{d}x}\approx \frac{h}{3}\left(\begin{array}{l}1\underset{\left(n-1\right)\text{-componentes}}{\underbrace{42424\cdots 4}}1\hfill \end{array}\right)\left(\begin{array}{c}f\left(x_0\right)\\ f\left(x_1\right)\\ f\left(x_2\right)\\ f\left(x_3\right)\\ f\left(x_4\right)\\ ⋮\\ f\left(x_n\right)\end{array}\right)$ (4)

To ensure the applicability of the composite rule, we assume that $n=2k$ , with $k\ge 2$ . The smallest composite case occurs when $k=2$ ; for $k=1$ , the formula reduces to the simple Simpson’s 1/3 rule.

2.3. Simpson’s 1/3 Rule and the Matrix Formulation for Composite Double Integration over a Rectangular Domain with Uniformly Spaced Data

We now derive the matrix representation of the double integral by applying the composite Simpson’s 1/3 rule twice, as defined in Equation (2):

$\begin{array}{l}{\displaystyle {\int }_c^d{\displaystyle {\int }_a^{b}f\left(x,y\right)\text{d}x\text{d}y}}\\ \approx {\displaystyle {\int }_c^d\frac{h_1}{3}}\left[f\left(x_0,y\right)+4{\displaystyle \sum _{\begin{array}{c}i=1\\ \text{Δ}i=2\end{array}}^{n-1}f}\left(x_i,y\right)+2{\displaystyle \sum _{\begin{array}{c}i=2\\ \text{Δ}i=2\end{array}}^{n-2}f}\left(x_i,y\right)+f\left(x_n,y\right)\right]\text{d}y\end{array}$

Let

$g\left(y\right)=f\left(x_0,y\right)+4{\displaystyle \sum _{\begin{array}{c}i=1\\ \text{Δ}i=2\end{array}}^{n-1}f}\left(x_i,y\right)+2{\displaystyle \sum _{\begin{array}{c}i=2\\ \text{Δ}i=2\end{array}}^{n-2}f}\left(x_i,y\right)+f\left(x_n,y\right)$

Then

$\begin{array}{c}{\displaystyle {\int }_c^d{\displaystyle {\int }_a^{b}f\left(x,y\right)\text{d}x\text{d}y}}\approx \frac{h_1}{3}{\displaystyle {\int }_c^{d}g\left(y\right)\text{d}y}\\ \approx \frac{h_1{h}_2}{9}\left[g\left(y_0\right)+4{\displaystyle \sum _{\begin{array}{c}i=1\\ \text{Δ}i=2\end{array}}^{m-1}g}\left(y_i\right)+2{\displaystyle \sum _{\begin{array}{c}i=2\\ \text{Δ}i=2\end{array}}^{m-2}g}\left(y_i\right)+g\left(y_m\right)\right]\end{array}$ (5)

Here, $h_1=\frac{b-a}{n}$ and $h_2=\frac{d-c}{m}$ are the step sizes along the $x$ - and $y$ - axes, respectively, and

$g\left(y_j\right)=f\left(x_0,y_j\right)+4{\displaystyle \sum _{\begin{array}{c}i=1\\ \text{Δ}i=2\end{array}}^{n-1}f}\left(x_i,y_j\right)+2{\displaystyle \sum _{\begin{array}{c}i=2\\ \text{Δ}i=2\end{array}}^{n-2}f}\left(x_i,y_j\right)+f\left(x_n,y_j\right)$

This expression can be written in matrix form as

$g\left(y_j\right)=\left(\begin{array}{l}1424\cdots 241\hfill \end{array}\right)\left(\begin{array}{c}f\left(x_0,y_j\right)\\ f\left(x_1,y_j\right)\\ f\left(x_2,y_j\right)\\ f\left(x_3,y_j\right)\\ f\left(x_4,y_j\right)\\ ⋮\\ f\left(x_n,y_j\right)\end{array}\right)$ (6)

On the other hand, Equation (5) can be expressed as

${\displaystyle {\int }_c^d{\displaystyle {\int }_a^{b}f\left(x,y\right)\text{d}x\text{d}y}}\approx \frac{h_1{h}_2}{9}\left(\begin{array}{l}g\left(y_0\right)g\left(y_1\right)g\left(y_2\right)\cdots g\left(y_m\right)\hfill \end{array}\right)\left(\begin{array}{c}1\\ 4\\ 2\\ 4\\ 2\\ 4\\ ⋮\\ 4\\ 1\end{array}\right)$ (7)

Taking Equation (6) into account, Equation (7) can be written compactly as

${\displaystyle {\int }_c^d{\displaystyle {\int }_a^{b}f\left(x,y\right)\text{d}x\text{d}y}}\approx \frac{h_1{h}_2}{9}KML$ (8)

where

$K={\left(\begin{array}{llllllll}1\hfill & 4\hfill & 2\hfill & 4\hfill & \cdots \hfill & 2\hfill & 4\hfill & 1\hfill \end{array}\right)}_{1\times \left(n+1\right)},L={\left(\begin{array}{c}1\\ 4\\ 2\\ 4\\ ⋮\\ 2\\ 4\\ 1\end{array}\right)}_{\left(m+1\right)\times 1}$

and

$M=\left(\begin{array}{ccccc}f\left(x_0,y_0\right)& f\left(x_0,y_1\right)& f\left(x_0,y_2\right)& \cdots & f\left(x_0,y_m\right)\\ f\left(x_1,y_0\right)& f\left(x_1,y_1\right)& f\left(x_1,y_2\right)& \cdots & f\left(x_1,y_m\right)\\ f\left(x_2,y_0\right)& f\left(x_2,y_1\right)& f\left(x_2,y_2\right)& \cdots & f\left(x_2,y_m\right)\\ ⋮& ⋮& ⋮& & ⋮\\ f\left(x_n,y_0\right)& f\left(x_n,y_1\right)& f\left(x_n,y_2\right)& \cdots & f\left(x_n,y_m\right)\end{array}\right)$ (9)

For composite integration, we require $n=2l$ , $l\ge 2$ and $m=2k$ , $k\ge 2$ . The minimal composite case corresponds to $l=2$ and $k=2$ . When $l=1$ and $k=1$ , the formula reduces to the simple (non-composite) Simpson’s rule for double integration.

The matrix M in Equation (9), is associated with the function values over the rectangular grid defining the double integral. Up to this point, the formulation is classical and well established in the literature. It is evident that, as the number of nodes increases, the size of the matrix grows accordingly. In what follows, we develop a new approach showing that this matrix can always be reduced to a 3 × 3 matrix whose entries involve double summations.

Finally, if $x_i=a+i{h}_1$ , $i=0,1,2,\cdots ,n$ , $y_j=c+j{h}_2$ , $j=0,1,2,\cdots ,m$ , then the matrix M in Equation (9) takes the explicit form

$\left(\begin{array}{ccccc}f\left(a,c\right)& f\left(a,c+h_2\right)& f\left(a,c+2h_2\right)& \cdots & f\left(a,c+m{h}_2\right)\\ f\left(a+h_1,c\right)& f\left(a+h_1,c+h_2\right)& f\left(a+h_1,c+2h_2\right)& \cdots & f\left(a+h_1,c+m{h}_2\right)\\ f\left(a+2h_1,c\right)& f\left(a+2h_1,c+h_2\right)& f\left(a+2h_1,c+2h_2\right)& \cdots & f\left(a+2h_1,c+m{h}_2\right)\\ ⋮& ⋮& ⋮& & ⋮\\ f\left(a+n{h}_1,c\right)& f\left(a+n{h}_1,c+h_2\right)& f\left(a+2h_1,c+2h_2\right)& \cdots & f\left(a+n{h}_1,c+m{h}_2\right)\end{array}\right)$

(10)

Equation (8), together with M defined in Equation (10), represents the composite Simpson’s 1/3 rule for double integration in matrix form for uniformly spaced data.

2.4. Matrix Method for Non-Composite Double Integration Based on Simpson’s 1/3 Rule

To compute the double integral

${\displaystyle {\int }_c^d{\displaystyle {\int }_a^{b}f\left(x,y\right)\text{d}x\text{d}y}}$

using a non-composite approach based on Simpson’s 1/3 rule, we set $n=m=2$ in Equation (10). The associated matrix is then given by

$M=\left(\begin{array}{ccc}f\left(a,c\right)& f\left(a,c+h_2\right)& f\left(a,c+2h_2\right)\\ f\left(a+h_1,c\right)& f\left(a+h_1,c+h_2\right)& f\left(a+h_1,c+2h_2\right)\\ f\left(a+2h_1,c\right)& f\left(a+2h_1,c+h_2\right)& f\left(a+2h_1,c+2h_2\right)\end{array}\right)$ (11)

Substituting this matrix into Equation (8) yields the approximation of the double integral.

Figure 1 provides a geometric illustration of this particular case. The nodes are defined as

$x_0=a,x_1=a+h_1,x_2=a+2h_1,$

$y_0=c,y_1=c+h_2,y_2=c+2h_2,$

with step sizes $h_1=\frac{b-a}{2},h_2=\frac{d-c}{2}$ . As an application, consider the double integral

Figure 1. (a) Nodes $x_0,x_1,x_2$ and $y_0,y_1,y_2$ define the rectangular domain $\left[x_0,x_2\right]\times \left[y_0,y_2\right]$ . (b) The surface $z=f\left(x,y\right)$ (blue) is approximated by an interpolated surface (red) with parabolic profiles through the nodes. The blue section, parallel to the $xz$ -plane, represents Simpson’s 1/3 integration along the $x$ -axis, while the red arrow indicates integration along the $y$ -axis with step size $h_2$ , yielding a volume for $f\left(x,y\right)\ge 0$ . This illustrates the geometric interpretation of the double integral approximation.

${\displaystyle {\int }_1^3{\displaystyle {\int }_0^4{\text{e}}^{x+y}\text{d}x\text{d}y}}$

In this case, $a=0$ , $b=4$ ; $c=1$ , $d=3$ . The step sizes are

$h_1=\frac{b-a}{2}=\frac{4-0}{2}=2;h_2=\frac{d-c}{2}=\frac{3-1}{2}=1$

and $f\left(x,y\right)={\text{e}}^{x+y}$ .

From Equation (11), the corresponding matrix is

$M=\left(\begin{array}{ccc}{\text{e}}^{a+c}& {\text{e}}^{a+\left(c+h_2\right)}& {\text{e}}^{a+\left(c+2h_2\right)}\\ {\text{e}}^{\left(a+h_1\right)+c}& {\text{e}}^{\left(a+h_1\right)+\left(c+h_2\right)}& {\text{e}}^{\left(a+h_1\right)+\left(c+2h_2\right)}\\ {\text{e}}^{\left(a+2h_1\right)+c}& {\text{e}}^{\left(a+2h_1\right)+\left(c+h_2\right)}& {\text{e}}^{\left(a+2h_1\right)+\left(c+2h_2\right)}\end{array}\right)$

Substituting $a=0$ , $c=1$ , $h_1=2$ and $h_2=1$ , we obtain

$M=\left(\begin{array}{ccc}{\text{e}}^1& {\text{e}}^2& {\text{e}}^3\\ {\text{e}}^3& {\text{e}}^4& {\text{e}}^5\\ {\text{e}}^5& {\text{e}}^6& {\text{e}}^7\end{array}\right)$

Therefore, the Simpson’s 1/3 approximation for the double integral is

${\displaystyle {\int }_1^3{\displaystyle {\int }_0^4{\text{e}}^{x+y}\text{d}x\text{d}y}}\approx \left(\frac{\left(2\right)\left(1\right)}{9}\right)\left(\begin{array}{lll}1\hfill & 4\hfill & 1\hfill \end{array}\right)\left(\begin{array}{ccc}{\text{e}}^1& {\text{e}}^2& {\text{e}}^3\\ {\text{e}}^3& {\text{e}}^4& {\text{e}}^5\\ {\text{e}}^5& {\text{e}}^6& {\text{e}}^7\end{array}\right)\left(\begin{array}{c}1\\ 4\\ 1\end{array}\right)\approx 990.82$

The exact value of the integral is 930.853; hence, the relative percentage error of the approximation is 6.44.

We now derive, in general form, a matrix method for the double integral

${\displaystyle {\int }_c^d{\displaystyle {\int }_a^{b}f\left(x,y\right)\text{d}x\text{d}y}}$

using Simpson’s 1/3 matrix rule when the nodes are not uniformly spaced. Along the $x$ -axis, the nodes are $a=x_0,x_2,x_4,\cdots ,x_n=b$ , and along the $y$ -axis, the nodes are $c=y_0,y_2,y_4,\cdots ,y_m=d$ . For convenience, we consider only nodes with even subscripts. The nodes with odd subscripts arise as midpoints between each pair of consecutive even nodes, thereby allowing the application of Simpson’s 1/3 rule. This construction generates a set of subrectangles defined by consecutive pairs of even nodes along both the $x$ - and $y$ - axes. Within each subrectangle, Equation (11) is employed to apply the matrix formulation of Simpson´s 1/3 rule for double integration. By summing over all such subrectangles, we obtain a composite integration formula suitable for non-uniformly spaced data, as will be shown in the next section.

In what follows, we build on the development presented in Section 2.4 to derive the general matrix formulation for composite numerical double integration.

3. A 3 × 3 Matrix-Based Method for General Composite Double Integration with Non-Uniformly Spaced Data Using Simpson’s 1/3 Rule

3.1. Formulation of the 3 × 3 Matrix-Based Method

To numerically compute the double integral

${\displaystyle {\int }_c^d{\displaystyle {\int }_a^{b}f\left(x,y\right)\text{d}x\text{d}y}}$

we apply the matrix formulation of the Simpson’s 1/3 rule for double integration, with the associated matrix M given in Equation (11). To this end, the domain $\left[a,b\right]\times \left[c,d\right]$ is partitioned into subrectangles of the form

$\left[x_i,x_{i+2}\right]\times \left[y_j,y_{j+2}\right]\subset \left[a,b\right]\times \left[c,d\right]$

where $i$ and $j$ are even indices, i.e., $i=0,2,\cdots ,n-2$ and $j=0,2,\cdots ,m-2$ .

The nodes along the x- and y-axes are arbitrarily prescribed and non-uniformly spaced. Since Simpson’s 1/3 rule is used, midpoint data are required within each subrectangle in both coordinate directions. These midpoints are defined as

$x_{i+1}=x_i+h_1^{\frac{i+2}{2}},y_{j+1}=y_j+h_1^{\frac{j+2}{2}}$

where $h_1^{\frac{i+2}{2}}$ and $h_2^{\frac{j+2}{2}}$ denote the local half-step sizes, given by

$h_1^{\frac{i+2}{2}}=\frac{x_{i+2}-x_i}{2},h_2^{\frac{j+2}{2}}=\frac{x_{j+2}-x_j}{2}$ (12)

(See Figure 2).

Figure 2. The nodes $x_i$ , with $i=2r$ , $r=0,\cdots ,\frac{n}{2}$ , along the $x$ -axis, and the nodes $y_j$ , with $j=2s$ , $s=0,\cdots ,\frac{m}{2}$ , along the $y$ - axis, are not uniformly spaced and define a rectangular mesh. The nodes with odd subscripts, $x_{i+1}$ and $y_{j+1}$ , correspond to the midpoints of the intervals $\left[x_i,x_{i+2}\right]$ and $\left[y_j,y_{j+2}\right]$ , respectively.

Let $I_{i\left(i+2\right);j\left(j+2\right)}$ denote the numerical approximation of the integral over the subrectangle $\left[x_i,x_{i+2}\right]\times \left[y_j,y_{j+2}\right]$ . By applying the matrix form of Simpson’s 1/3 rule to each subrectangle, we obain

$I_{i\left(i+2\right);j\left(j+2\right)}=\frac{\left(x_{i+2}-x_i\right)\left(y_{j+2}-y_j\right)}{36}\left(\begin{array}{lll}1\hfill & 4\hfill & 1\hfill \end{array}\right)M_{ij}\left(\begin{array}{c}1\\ 4\\ 1\end{array}\right)$ (13)

where

$M_{ij}=\left(\begin{array}{ccc}f\left(x_i,y_j\right)& f\left(x_i,y_j+h_2^{\frac{j+2}{2}}\right)& f\left(x_i,y_j+2h_2^{\frac{j+2}{2}}\right)\\ f\left(x_i+h_1^{\frac{i+2}{2}},y_j\right)& f\left(x_i+h_1^{\frac{i+2}{2}},y_j+h_2^{\frac{j+2}{2}}\right)& f\left(x_i+h_1^{\frac{i+2}{2}},y_j+2h_2^{\frac{j+2}{2}}\right)\\ f\left(x_i+2h_1^{\frac{i+2}{2}},y_j\right)& f\left(x_i+2h_1^{\frac{i+2}{2}},y_j+h_2^{\frac{j+2}{2}}\right)& f\left(x_i+2h_1^{\frac{i+2}{2}},y_j+2h_2^{\frac{j+2}{2}}\right)\end{array}\right)$

Which is consistent with Equation (11) in Section 2.4.

The composite approximation to the double integral over $\left[a,b\right]\times \left[c,d\right]$ is obtained by summing over all subrectangles:

${\displaystyle {\int }_c^d{\displaystyle {\int }_a^{b}f\left(x,y\right)\text{d}x\text{d}y}}\approx {\displaystyle \sum _{\begin{array}{c}i=0\\ \text{Δ}i=2\end{array}}^{n-2}{\displaystyle \sum _{\begin{array}{c}j=0\\ \text{Δ}j=2\end{array}}^{m-2}{I}_{i\left(i+2\right);j\left(j+2\right)}}}$

To derive the global matrix form, we first consider representative terms for different values of $i$ and $j$ in $I_{i\left(i+2\right);j\left(j+2\right)}$ , for example:

$I_{02;02},I_{02;24},I_{02;46},I_{24;02},I_{24;24},I_{24;46},I_{46;02},I_{46;24},I_{46;46}$

each having the structure given in Equation (13) with the corresponding indices. Thus, the approximation to the integral is given by the sum of all such terms:

$\begin{array}{c}{\displaystyle {\int }_c^d{\displaystyle {\int }_a^{b}f\left(x,y\right)\text{d}x\text{d}y}}\approx I_{02;02}+I_{02;24}+I_{02;46}+I_{24;02}+I_{24;24}+I_{24;46}\\ +I_{46;02}+I_{46;24}+I_{46;46}+\cdots \end{array}$

This expression can be written compactly in matrix form as

${\displaystyle {\int }_c^d{\displaystyle {\int }_a^{b}f\left(x,y\right)\text{d}x\text{d}y}}\approx \left(\begin{array}{lll}1\hfill & 4\hfill & 1\hfill \end{array}\right)M\left(\begin{array}{c}1\\ 4\\ 1\end{array}\right)$ (14)

where $M$ is a 3 × 3 matrix whose elements are defined as follows

$M_{11}={\displaystyle \sum _{\begin{array}{c}j=0\\ \text{Δ}j=2\end{array}}^{m-2}{\displaystyle \sum _{\begin{array}{c}i=0\\ \text{Δ}i=2\end{array}}^{n-2}\frac{\left(x_{i+2}-x_i\right)\left(y_{j+2}-y_j\right)}{36}f\left(x_i,y_j\right)}}$

$M_{21}={\displaystyle \sum _{\begin{array}{c}j=0\\ \text{Δ}j=2\end{array}}^{m-2}{\displaystyle \sum _{\begin{array}{c}i=0\\ \text{Δ}i=2\end{array}}^{n-2}\frac{\left(x_{i+2}-x_i\right)\left(y_{j+2}-y_j\right)}{36}f\left(x_i+h_1^{\frac{i+2}{2}},y_j\right)}}$

$M_{31}={\displaystyle \sum _{\begin{array}{c}j=0\\ \text{Δ}j=2\end{array}}^{m-2}{\displaystyle \sum _{\begin{array}{c}i=0\\ \text{Δ}i=2\end{array}}^{n-2}\frac{\left(x_{i+2}-x_i\right)\left(y_{j+2}-y_j\right)}{36}f\left(x_i+2h_1^{\frac{i+2}{2}},y_j\right)}}$

$M_{12}={\displaystyle \sum _{\begin{array}{c}j=0\\ \text{Δ}j=2\end{array}}^{m-2}{\displaystyle \sum _{\begin{array}{c}i=0\\ \text{Δ}i=2\end{array}}^{n-2}\frac{\left(x_{i+2}-x_i\right)\left(y_{j+2}-y_j\right)}{36}f\left(x_i,y_j+h_2^{\frac{j+2}{2}}\right)}}$

$M_{22}={\displaystyle \sum _{\begin{array}{c}j=0\\ \text{Δ}j=2\end{array}}^{m-2}{\displaystyle \sum _{\begin{array}{c}i=0\\ \text{Δ}i=2\end{array}}^{n-2}\frac{\left(x_{i+2}-x_i\right)\left(y_{j+2}-y_j\right)}{36}f\left(x_i+h_1^{\frac{i+2}{2}},y_j+h_2^{\frac{j+2}{2}}\right)}}$ (15)

$M_{32}={\displaystyle \sum _{\begin{array}{c}j=0\\ \text{Δ}j=2\end{array}}^{m-2}{\displaystyle \sum _{\begin{array}{c}i=0\\ \text{Δ}i=2\end{array}}^{n-2}\frac{\left(x_{i+2}-x_i\right)\left(y_{j+2}-y_j\right)}{36}f\left(x_i+2h_1^{\frac{i+2}{2}},y_j+h_2^{\frac{j+2}{2}}\right)}}$

$M_{13}={\displaystyle \sum _{\begin{array}{c}j=0\\ \text{Δ}j=2\end{array}}^{m-2}{\displaystyle \sum _{\begin{array}{c}i=0\\ \text{Δ}i=2\end{array}}^{n-2}\frac{\left(x_{i+2}-x_i\right)\left(y_{j+2}-y_j\right)}{36}f\left(x_i,y_j+2h_2^{\frac{j+2}{2}}\right)}}$

$M_{23}={\displaystyle \sum _{\begin{array}{c}j=0\\ \text{Δ}j=2\end{array}}^{m-2}{\displaystyle \sum _{\begin{array}{c}i=0\\ \text{Δ}i=2\end{array}}^{n-2}\frac{\left(x_{i+2}-x_i\right)\left(y_{j+2}-y_j\right)}{36}f\left(x_i+h_1^{\frac{i+2}{2}},y_j+2h_2^{\frac{j+2}{2}}\right)}}$

$M_{33}={\displaystyle \sum _{\begin{array}{c}j=0\\ \text{Δ}j=2\end{array}}^{m-2}{\displaystyle \sum _{\begin{array}{c}i=0\\ \text{Δ}i=2\end{array}}^{n-2}\frac{\left(x_{i+2}-x_i\right)\left(y_{j+2}-y_j\right)}{36}f\left(x_i+2h_1^{\frac{i+2}{2}},y_j+2h_2^{\frac{j+2}{2}}\right)}}$

Equation (14) and Equation (15) together constitute the formulation of the 3 × 3 matrix-based method.

Equation (15) can be further simplified by assuming uniform spacing; that is,

$h=h_1^k=h_2^l,k,l=1,2,3,\cdots ,\frac{n}{2},\frac{m}{2}$

and

$H=x_{i+2}-x_i=y_{j+2}-y_j,i,j=0,1,2,\cdots ,\left(n-2\right),\left(m-2\right)$

Under these assumptions, the expressions above simplify considerably due to the constant mesh width.

For a numerical equivalence of the matrix-based relation to verify the validity of the relation

${\displaystyle \sum _{\begin{array}{c}i=0\\ \text{Δ}i=2\end{array}}^{n-2}{\displaystyle \sum _{\begin{array}{c}j=0\\ \text{Δ}j=2\end{array}}^{m-2}{I}_{i\left(i+2\right);j\left(j+2\right)}}}=\left(\begin{array}{lll}1\hfill & 4\hfill & 1\hfill \end{array}\right)M\left(\begin{array}{c}1\\ 4\\ 1\end{array}\right)$

we consider a specific numerical example. The purpose of this example is to show that the sum of the local composite Simpson contributions over all subrectangles is equal to the corresponding matrix-based expression.

Consider the double integral

${\displaystyle {\int }_0^2{\displaystyle {\int }_0^1{\text{e}}^{x+y}\text{d}x\text{d}y}}$

For the numerical evaluation, the prescribed nodes are

For the $x$ -direction For the $y$ -direction

$x_0=0,x_2=0.4,x_4=1$ $y_0=0,y_2=1.2,y_4=2$

Thus, $n=m=4$ and according to Equation (12), $h_1^1=0.2$ , $h_1^2=0.3$ , $h_2^1=0.6$ , $h_2^2=0.4$ . The corresponding partition contains four local Simpson subrectangles. Therefore, the double sum can be expanded as

$I_{02;02},I_{02;24},I_{24;02},I_{24;24}$

We first evaluate each local contribution separately using Equation (13). For the first subrectangle, corresponding to $I_{02;02}$ , we obtain

$I_{02;02}=\frac{\left(x_2-x_0\right)\left(y_2-y_0\right)}{36}\left(\begin{array}{lll}1\hfill & 4\hfill & 1\hfill \end{array}\right)\left(\begin{array}{ccc}f\left(x_0,y_0\right)& f\left(x_0,y_0+h_2^1\right)& f\left(x_0,y_0+2h_2^1\right)\\ f\left(x_0+h_1^1,y_0\right)& f\left(x_0+h_1^1,y_0+h_2^1\right)& f\left(x_0+h_1^1,y_0+2h_2^1\right)\\ f\left(x_0+2h_1^1,y_0\right)& f\left(x_0+2h_1^1,y_0+h_2^1\right)& f\left(x_0+2h_1^1,y_0+2h_2^1\right)\end{array}\right)\left(\begin{array}{c}1\\ 4\\ 1\end{array}\right)$

Substituting the prescribed nodes and $f\left(x,y\right)={\text{e}}^{x+y}$ , this gives

$I_{02;02}=\frac{\left(0.4\right)\left(1.2\right)}{36}\left(\begin{array}{lll}1\hfill & 4\hfill & 1\hfill \end{array}\right)\left(\begin{array}{ccc}f\left(0,0\right)& f\left(0,0.6\right)& f\left(0,1.2\right)\\ f\left(0.2,0\right)& f\left(0.2,0.6\right)& f\left(0.2,1.2\right)\\ f\left(0.4,0\right)& f\left(0.4,0.6\right)& f\left(0.4,1.2\right)\end{array}\right)\left(\begin{array}{l}1\hfill \\ 4\hfill \\ 1\hfill \end{array}\right)$

$I_{02;02}=\frac{\left(0.4\right)\left(1.2\right)}{36}\left(\begin{array}{lll}1\hfill & 4\hfill & 1\hfill \end{array}\right)\left(\begin{array}{ccc}1& {\text{e}}^{0.6}& {\text{e}}^{1.2}\\ {\text{e}}^{0.2}& {\text{e}}^{0.8}& {\text{e}}^{1.4}\\ {\text{e}}^{0.4}& {\text{e}}^1& {\text{e}}^{1.6}\end{array}\right)\left(\begin{array}{l}1\hfill \\ 4\hfill \\ 1\hfill \end{array}\right)$

$I_{02;02}=1.14188856421427$

Similarly, for the rest of the subrectangles $I_{02;24}$ , $I_{24;02}$ ,$I_{24;24}$

$I_{02;24}=2.00150177080648$

$I_{24;02}=2.84761489484506$

$I_{24;24}=4.99129812945376$

Adding the four local contributions yields

${\displaystyle \sum _{\begin{array}{c}i=0\\ \text{Δ}i=2\end{array}}^{n-2}{\displaystyle \sum _{\begin{array}{c}j=0\\ \text{Δ}j=2\end{array}}^{m-2}{I}_{i\left(i+2\right);j\left(j+2\right)}}}=10.9823033593196$

We now verify the same result using the matrix representation. The entries of the matrix M are obtained by collecting the corresponding function values from the four local Simpson subrectangles. For example,

$\begin{array}{c}{M}_{11}=\frac{\left(x_2-x_0\right)\left(y_2-y_0\right)}{36}f\left(x_0,y_0\right)+\frac{\left(x_4-x_2\right)\left(y_2-y_0\right)}{36}f\left(x_2,y_0\right)\\ +\frac{\left(x_2-x_0\right)\left(y_4-y_2\right)}{36}f\left(x_0,y_2\right)+\frac{\left(x_4-x_2\right)\left(y_4-y_2\right)}{36}f\left(x_2,y_2\right)\end{array}$

Substitution of the prescribed nodes and $f\left(x,y\right)={\text{e}}^{x+y}$ , gives

$\begin{array}{c}{M}_{11}=\frac{\left(0.4\right)\left(1.2\right)}{36}f\left(0,0\right)+\frac{\left(0.6\right)\left(1.2\right)}{36}f\left(0.4,0\right)\\ +\frac{\left(0.4\right)\left(0.8\right)}{36}f\left(0,1.2\right)+\frac{\left(0.6\right)\left(0.8\right)}{36}f\left(0.4,1.2\right)\end{array}$

$M_{11}=\frac{\left(0.4\right)\left(1.2\right)}{36}+\frac{\left(0.6\right)\left(1.2\right)}{36}{\text{e}}^{0.4}+\frac{\left(0.4\right)\left(0.8\right)}{36}{\text{e}}^{1.2}+\frac{\left(0.6\right)\left(0.8\right)}{36}{\text{e}}^{1.6}$

$M_{11}=0.138722410035752$

The remaining entries of $M$ are obtained in the same manner. Consequently,

$M=\left(\begin{array}{ccc}0.138722410035752& 0.221208256776975& 0.355985057856772\\ 0.181751905416110& 0.289823555924596& 0.466405986952344\\ 0.238617100747032& 0.380501411991207& 0.612331651339934\end{array}\right)$

Therefore

$\left(\begin{array}{lll}1\hfill & 4\hfill & 1\hfill \end{array}\right)M\left(\begin{array}{c}1\\ 4\\ 1\end{array}\right)=10.9823033593196$

This value coincides with the value obtained by directly adding the four local Simpson contributions. Hence,

${\displaystyle \sum _{\begin{array}{c}i=0\\ \text{Δ}i=2\end{array}}^2{\displaystyle \sum _{\begin{array}{c}j=0\\ \text{Δ}j=2\end{array}}^2{I}_{i\left(i+2\right);j\left(j+2\right)}}}=\left(\begin{array}{lll}1\hfill & 4\hfill & 1\hfill \end{array}\right)M\left(\begin{array}{c}1\\ 4\\ 1\end{array}\right)$

which numerically verifies the stated relation for the considered non-uniformly spaced nodes.

3.2. Error Analysis of the 3 × 3 Matrix-Based Method

The integration error over the interval $\left[x_{k-1},x_{k+1}\right]$ using Simpson’s 1/3 is given (See Yang et al. [11]) by

$\begin{array}{l}{\displaystyle {\int }_{x_{k-1}}^{x_{k+1}}f\left(x\right)\text{d}x}-\frac{h}{3}\left(f\left(x_{k-1}\right)+4f\left(x_k\right)+f\left(x_{k+1}\right)\right)\\ =-\frac{1}{90}{h}^5{f}^{\left(4\right)}\left(x_k\right)+O\left(h^7\right)=O\left(h^5\right)\end{array}$ (16)

This shows that the local truncation error over $\left[x_{k-1},x_{k+1}\right]$ is proportional to $h^5$ . Therefore, the error over the entire interval $\left[a,b\right]$ , where $a=x_0$ and $b=x_n$ is,

${\displaystyle {\int }_{x_0}^{x_n}f\left(x\right)\text{d}x}-\frac{h}{3}{\displaystyle \sum }_{\begin{array}{c}k=1\\ \text{Δ}k=2\end{array}}^{n-1}\left(f\left(x_{k-1}\right)+4f\left(x_k\right)+f\left(x_{k+1}\right)\right)=\frac{n}{2}O\left(h^5\right)$ (17)