A Matrix-Based Reformulation of the Composite Simpson Rule for Double Numerical Integration over Rectangular Domains with Non-Uniformly Spaced Data

Lizandro B. R. Zegarra^1*, Luis E. G. Armas^2*, César C. Campos³, Pedro G. R. Carrera¹, Herón J. M. Marchena¹, Carlos A. M. Alvarado¹, María E. B. Gemin¹, Aída B. D. Reyna¹

¹Department of Mathematics, Universidad Nacional del Santa (UNS), Nuevo Chimbote, Peru.
²Optics and Micro/Nanofabrication of Devices Group (GOMNDI), Federal University of Pampa, Campus Alegrete, Rio Grande do Sul (RS), Brazil.
³Department of Mathematics, Universidad Nacional de Huancavelica (UNH), Huancavelica, Peru.

DOI: 10.4236/jamp.2026.1410200                

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

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)

Since $n$ is an even integer and $n~\frac{1}{h}$ , the global error is $O\left(h^4\right)$ .

For composite double integration over the rectangle $\left[x_0,x_n\right]\times \left[y_0,y_m\right]$ , there are $\frac{n}{2}$ subintervals along the $x$ -axis and $\frac{m}{2}$ subintervals along the $y$ -axis (see Figure 2), where $n$ , $m$ are even positive integers. Thus, Simpson’s 1/3 rule is applied $\frac{n}{2}\times \frac{m}{2}$ times. Consider the double integral over a subrectangle $\left[x_i,x_j\right]\times \left[y_k,y_l\right]$ , with $j-i=2$ and $l-k=2$

$I={\displaystyle {\int }_{y_k}^{y_l}{\displaystyle {\int }_{x_i}^{x_j}f\left(x,y\right)\text{d}x\text{d}y}}$

Let

$I_1\left(y\right)={\displaystyle {\int }_{x_i}^{x_j}f\left(x,y\right)\text{d}x}$

Applying Equation (16) with respect to $x$ , we obtain

$I_1\left(y\right)=\frac{h_1}{3}\left[f\left(x_i,y\right)+4f\left(x_{i+1},y\right)+f\left(x_{i+2},y\right)\right]-\frac{1}{90}{h}_1^5{f}_x^{\left(4\right)}\left(x_{i+1},y\right)+O\left(h_1^7\right)$

where $h_1$ is actually $h_1^{\frac{i+2}{2}}$ .

Applying Simpson’s rule again with respect to $y$ , we obtain

$I={\displaystyle {\int }_{y_k}^{y_l}{I}_1\left(y\right)\text{d}y}=\frac{h_2}{3}\left[I_1\left(y_k\right)+4I_1\left(y_{k+1}\right)+I_1\left(y_{k+2}\right)\right]-\frac{1}{90}{h}_2^5{I}_1^{\left(4\right)}\left(y_{k+1}\right)+O\left(h_2^7\right)$

where $h_2$ is actually $h_2^{\frac{i+2}{2}}$ .

Substituting $I_1\left(y\right)$ into this expression, shows that the leading truncation error over the subrectangle $\left[x_i,x_j\right]\times \left[y_k,y_l\right]$ , is proportional to $h_1{h}_2^5$ or $h_2{h}_1^5$ .

If $h_1$ and $h_2$ are of the same order (i.e., $h_1~h_2~h$ ) then the local integration error becomes $O\left(h^6\right)$ . Thus,

${\displaystyle {\int }_{y_k}^{y_l}{\displaystyle {\int }_{x_i}^{x_j}f\left(x,y\right)\text{d}x\text{d}y}}-\frac{h}{3}\left[I_{N_1}+4I{N}_2+I_{N_3}\right]=O\left(h^6\right)$

where

• $I_{N_1}$ is the Simpson approximation of ${\displaystyle {\int }_{y_k}^{y_l}f\left(x_i,y\right)\text{d}y}$

• $I_{N_2}$ is the Simpson approximation of ${\displaystyle {\int }_{y_k}^{y_l}f\left(x_{i+1},y\right)\text{d}y}$

• $I_{N_3}$ is the Simpson approximation of ${\displaystyle {\int }_{y_k}^{y_l}f\left(x_{i+2},y\right)\text{d}y}$

This shows that the local integration error over each subrectangle is $O\left(h^6\right)$ . Therefore, using Equation (17), the global integration error over the entire domain $\left[x_0,x_n\right]\times \left[y_0,y_m\right]$ is

${\displaystyle {\int }_{y_0}^{y_m}{\displaystyle {\int }_{x_0}^{x_n}f\left(x,y\right)\text{d}x\text{d}y}}-\left(\begin{array}{lll}1\hfill & 4\hfill & 1\hfill \end{array}\right)M\left(\begin{array}{c}1\\ 4\\ 1\end{array}\right)=\frac{nm}{4}O\left(h^6\right)=O\left(h^4\right)$

Hence, the global error of the method is of order $O\left(h^4\right)$ , since both factors $\frac{n}{2}$ and $\frac{m}{2}$ are proportional to $\frac{1}{h}$ .

3.3. Comparison with Gaussian Quadrature and Adaptive Cubature Methods

Gaussian quadrature and adaptive numerical integration methods are widely used because of their accuracy and flexibility, and both approaches can be extended from one-dimensional to multidimensional integrals. For double integrals over rectangular domains, one-dimensional Gaussian rules can be combined through tensor products, while adaptive strategies recursively subdivide the integration domain according to local error estimates.

Gaussian quadrature is particularly effective for analytically defined and sufficiently smooth functions, since its nodes and weights are determined by the selected quadrature rule. However, when data are prescribed at non-uniformly spaced locations, the Gaussian nodes generally do not coincide with the available data, and additional function evaluations or interpolation may be required. The proposed method is intended for this situation. It operates directly on prescribed non-uniform meshes without redistributing the available nodes and applies Simpson’s 1/3 rule locally over rectangular subdomains determined by consecutive pairs of nodes.

Adaptive multidimensional integration methods, commonly referred to as adaptive cubature methods, estimate local integration errors and refine the mesh in regions where the integrand requires additional resolution. This strategy can provide high accuracy and automatic error control, although it involves additional procedures for error estimation, domain subdivision, and refinement.

The proposed method does not perform automatic mesh refinement or local error control. Instead, it assumes a prescribed non-uniform mesh and applies the composite Simpson’s 1/3 formulation through a fixed 3 × 3 matrix representation. Its main advantage is the compact organization of the integration through a fixed 3 × 3 matrix. Although the nine matrix entries are obtained from double summations over the subrectangles, the matrix size remains independent of the number of discretization nodes. This provides a compact representation and can reduce the auxiliary memory required by the matrix formulation.

Under the uniform-spacing assumptions discussed in Section 3.2, the proposed method retains the fourth-order behavior of the composite Simpson’s 1/3 rule. Therefore, the proposed formulation should be regarded as a complementary alternative to Gaussian quadrature and adaptive cubature, particularly when function values are available on a prescribed non-uniform rectangular mesh and a simple fixed-size matrix representation is desirable. Table 1 presents a summary comparison of numerical integration methods.

Table 1. Comparison of numerical integration methods.

Method	Node distribution	Adaptability	Typical accuracy	Main characteristic
Classical composite Simpson	Uniform	Low	Fourth order	Simple and effective for Smooth functions
Proposed method	Prescribed non-uniform mesh	Low-moderate	Fourth order for uniform spacing; convergence order not established for general non-uniform meshes	Fixed 3 × 3 matrix representation
Gaussian quadrature	Nodes selected by the rule	low	High for smooth functions	Efficient use of optimally selected nodes
Adaptive cubature	Dynamically selected refined	High	Depends on tolerance and method	Automatic refinement based on error estimates

Throughout the numerical examples, the term prescribed nodes refer to the nodes with even indices in each coordinate direction. The intervening odd-indexed nodes are midpoint nodes determined by Equation (12). Thus, if N prescribed nodes are used in each coordinate direction, the resulting mesh contains 2N − 1 nodes in each direction, including the midpoint nodes.

4. Examples

As a first example, we numerically evaluate the double integral

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

which was previously computed using the method presented in Section 2.4. In this case, we consider non-uniformly spaced nodes prescribed as

$a=x_0=0,x_2=1.4,x_4=2.9,b=x_6=4,$

$c=y_0=1,y_2=1.9,y_4=2.6,d=y_6=3$

Thus, $n=6$ and $m=6$ . Using Equation (12), the corresponding local half-step sizes are computed as follows:

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

$\begin{array}{l}i=0:h_1^1=\frac{x_2-x_0}{2}=\frac{1.4-0}{2}=0.7\\ i=2:h_1^2=\frac{x_4-x_2}{2}=\frac{2.9-1.4}{2}=0.75\\ i=4:h_1^3=\frac{x_6-x_4}{2}=\frac{4-2.9}{2}=0.55\end{array}$ $\begin{array}{l}j=0:h_2^1=\frac{y_2-y_0}{2}=\frac{1.9-1}{2}=0.45\\ j=2:h_2^2=\frac{y_4-y_2}{2}=\frac{2.6-1.9}{2}=0.35\\ j=4:h_2^3=\frac{y_6-y_4}{2}=\frac{3-2.6}{2}=0.2\end{array}$

Next, the entries of the matrix M are computed according to Equation (15). Then, by applying Equation (14), we obtain

$I\approx \left(\begin{array}{lll}1\hfill & 4\hfill & 1\hfill \end{array}\right)\left(\begin{array}{ccc}9.54892362& 13.01679981& 17.91643942\\ 17.49847499& 23.85338444& 32.83201120\\ 32.32163104& 44.05985614& 60.64437918\end{array}\right)\left(\begin{array}{c}1\\ 4\\ 1\end{array}\right)\approx 931.71409288$

This result shows excellent agreement with the exact value, 930.85274423, yielding a relative percentage error of 0.09. Figure 3 depicts the function ${\text{e}}^{x+y}$ and Table 2 shows the convergence of the numerical integral. Using 40 prescribed nodes in each coordinate direction, the relative percentage error with respect to the exact value is 0.000014.

Figure 3. Graph of the function ${\text{e}}^{x+y}$ over $\left[0,4\right]\times \left[1,3\right]$ . The surface exhibits smooth, non-oscillatory behavior.

As another example, consider de evaluation of the double integral

${\displaystyle {\int }_0^3{\displaystyle {\int }_0^{3}xy\cos \left(x^2+y^2\right)\text{d}x\text{d}y}}$

using non-uniformly spaced nodes, prescribed as

$a=x_0=0,x_2=0.5,x_4=0.8,x_6=1.7,x_8=2.4,b=x_{10}=3$

$c=y_0=0,y_2=0.6,y_4=1.1,y_6=1.8,y_8=2.5,d=y_{10}=3$

Hence, $n=m=10$ . From Equation (12), the local half-step sizes are obtained as follows:

Table 2. Numerical convergence of the integral ${\displaystyle {\int }_1^3{\displaystyle {\int }_0^4{\text{e}}^{x+y}\text{d}x\text{d}y}}$ with exact value 930.85274423.

Number of prescribed nodes		Node values	Numerical value	Relative Error (%)
2	x:	0.0 4.0	990.81927095	6.44
	y:	1.0 3.0
4	x:	0.0 1.4 2.9 4.0	931.71409288	0.09
	y:	1.1 1.9 2.6 3.0
6	x:	0.0 0.5 0.9 1.8 2.7 4.0	931.55997772	0.08
	y	1.0 1.5 1.9 2.3 2.7 3.0
12	x:	0.0 0.4 0.7 1.0 1.4 1.7 2 2.5 2.8 3.0 3.5 4.0	930.86990692	0.002
	y:	1.0 1.1 1.3 1.5 1.9 2 2.2 2.5 2.6 2.8 2.9 3.0
16	x:	0.00 0.10 0.40 0.70 1.00 1.30 1.40 1.70 2.00 2.30 2.50 2.80 3.00 3.50 3.70 4.00	930.85937002	0.0007
	y:	1.00 1.15 1.35 1.45 1.60 1.80 2.00 2.20 2.30 2.50 2.57 2.60 2.80 2.86 2.90 3.00
20	x:	0.00 0.10 0.25 0.40 0.55 0.70 1.00 1.30 1.40 1.70 2.00 2.15 2.30 2.50 2.65 2.80 3.00 3.50 3.70 4.00	930.85889034	0.00066
	y:	1.00 1.15 1.25 1.35 1.45 1.60 1.75 1.80 1.95 2.00 2.20 2.30 2.45 2.50 2.57 2.60 2.80 2.86 2.90 3.00
40	x:	0.00 0.11 0.20 0.33 0.40 0.50 0.60 0.70 0.82 0.90 1.00 1.10 1.22 1.30 1.40 1.50 1.60 1.70 1.80 1.93 2.00 2.10 2.20 2.30 2.40 2.50 2.60 2.70 2.80 2.90 3.00 3.10 3.21 3.30 3.40 3.50 3.62 3.71 3.80 4.00	930.85287232	0.000014
	y:	1.00 1.05 1.10 1.15 1.22 1.25 1.30 1.35 1.40 1.45 1.51 1.55 1.60 1.65 1.70 1.75 1.80 1.85 1.93 1.95 2.02 2.05 2.12 2.15 2.20 2.25 2.31 2.35 2.40 2.45 2.53 2.55 2.60 2.65 2.70 2.75 2.83 2.88 2.96 3.00

For the $x$ -direction:

$h_1^1=0.25,h_1^2=0.15,h_1^3=0.45,h_1^4=0.35,h_1^5=0.30$

For the $y$ -direction:

$h_2^1=0.30,h_2^2=0.25,h_2^3=0.35,h_2^4=0.25,h_2^5=0.35$

Once the elements of the matrix $M$ , defined in Equation (15), have been computed, Equation (14) is applied to obtain:

$I\approx \left(\begin{array}{lll}1\hfill & 4\hfill & 1\hfill \end{array}\right)\left(\begin{array}{ccc}0.01164201& 0.00432489& -0.00904631\\ -0.01713386& -0.03333385& -0.04431376\\ -0.02771866& -0.02176820& -0.00297294\end{array}\right)\left(\begin{array}{c}1\\ 4\\ 1\end{array}\right)$

$I\approx -0.87700122$

Since the exact value is −0.87064431, the relative percentage error is 0.73.

Figure 4 shows $z=xy\cos \left(x^2+y^2\right)$ , and Table 3 reports the numerical convergence over $\left[0,3\right]\times \left[0,3\right]$ . Due to its oscillatory behavior, finer discretization is required. With 201 prescribed nodes, the relative percentage error with 5 decimal places is 0.00004.

Figure 4. Graph of $z=xy\cos \left(x^2+y^2\right)$ over $\left[0,3\right]\times \left[0,3\right]$ . Due to its oscillatory nature, a finer discretization is required to ensure convergence of the numerical integral.

Table 3. Numerical convergence of the integral ${\displaystyle {\int }_0^3{\displaystyle {\int }_0^{3}xy\cos \left(x^2+y^2\right)\text{d}x\text{d}y}}$ with exact value −0.87064431.

Number of prescribed nodes		Node values	Numerical value	Relative Error (%)
4	x:	0.0 1.2 2.0 3.0	−0.60395541	−30.63
	y:	0.0 1.0 2.2 3.0
16	x:	0.0 0.1 0.3 0.4 0.7 0.9 1.0 1.3 1.4 1.7 1.8 2.0 2.3 2.5 2.8 3.0	−0.87292387	0.26
	y:	0.0 0.1 0.2 0.4 0.7 0.9 1.0 1.3 1.4 1.7 1.8 2.1 2.3 2.5 2.7 3.0
20	x:	0.00 0.10 0.25 0.40 0.55 0.70 0.85 1.00 1.30 1.40 1.60 1.70 2.00 2.15 2.30 2.50 2.65 2.80 2.90 3.00	−0.87078921	0.017
	y:	0.00 0.10 0.25 0.35 0.55 0.70 0.85 1.00 1.30 1.40 1.60 1.85 2.00 2.15 2.30 2.50 2.60 2.75 2.90 3.00
40	x:	0.00 0.12 0.16 0.23 0.27 0.33 0.36 0.40 0.47 0.50 0.60 0.72 0.83 0.90 1.00 1.10 1.22 1.30 1.40 1.50 1.55 1.60 1.70 1.80 1.86 1.93 2.00 2.13 2.20 2.30 2.35 2.43 2.50 2.60 2.70 2.80 2.85 2.92 2.97 3.00	−0.87063454	−0.0011
	y:	0.00 0.11 0.15 0.20 0.27 0.33 0.36 0.40 0.47 0.50 0.60 0.70 0.82 0.90 1.00 1.10 1.22 1.30 1.40 1.50 1.55 1.60 1.70 1.80 1.86 1.93 2.00 2.10 2.20 2.30 2.35 2.40 2.50 2.60 2.70 2.80 2.85 2.90 2.96 3.00
201	x:	$x=\sqrt{z}$ , $0\le z\le 1$ , $z=0.045i$ , $i=0,1,\cdots ,200$	−0.87064462	0.00004
	y:	$y=\sqrt{z}$ , $0\le z\le 1$ , $z=0.045i$ , $i=0,1,\cdots ,200$

As a final example, consider the double integral

${\displaystyle {\int }_0^1{\displaystyle {\int }_0^1\sin \left(20x\right)+\sin \left(20y\right)\text{d}x\text{d}y}}$

Table 4 presents the numerical convergence over $\left[0,1\right]\times \left[0,1\right]$ . Due to the oscillatory integrand (Figure 5), finer discretization is required. With 801 prescribed nodes, the relative percentage error is 0.004.

Table 4. Convergence of the integral ${\displaystyle {\int }_0^1{\displaystyle {\int }_0^1\sin \left(20x\right)+\sin \left(20y\right)\text{d}x\text{d}y}}$ with exact value 0.05919179.

Number of prescribed nodes		Node values	Numerical value	Relative Error (%)
20	x:	0 0.10 0.15 0.19 0.25 0.30 0.34 0.41 0.45 0.56 0.60 0.65 0.69 0.73 0.79 0.85 0.90 0.95 0.97 1	0.05895511	−0.40
	y:	0 0.10 0.14 0.19 0.26 0.30 0.34 0.40 0.45 0.56 0.61 0.65 0.70 0.73 0.78 0.85 0.90 0.94 0.98 1
40	x:	0 0.10 0.12 0.13 0.15 0.18 0.21 0.23 0.25 0.28 0.30 0.34 0.37 0.38 0.41 0.44 0.48 0.50 0.56 0.58 0.60 0.62 0.65 0.66 0.69 0.72 0.73 0.75 0.77 0.79 0.81 0.82 0.85 0.88 0.90 0.92 0.95 0.97 0.99 1	0.05981197	1.05
	y:	0 0.09 0.12 0.14 0.15 0.19 0.21 0.23 0.25 0.28 0.30 0.34 0.37 0.38 0.41 0.44 0.48 0.50 0.56 0.58 0.60 0.62 0.65 0.66 0.70 0.71 0.73 0.76 0.77 0.79 0.81 0.82 0.84 0.88 0.90 0.92 0.95 0.97 0.98 1
60	x:	0.00 0.01 0.03 0.04 0.06 0.07 0.09 0.10 0.12 0.14 0.15 0.16 0.18 0.20 0.21 0.22 0.24 0.26 0.27 0.28 0.30 0.32 0.33 0.34 0.36 0.38 0.39 0.40 0.42 0.43 0.45 0.46 0.48 0.49 0.51 0.53 0.54 0.55 0.57 0.58 0.60 0.61 0.63 0.65 0.66 0.68 0.69 0.71 0.72 0.73 0.75 0.77 0.78 0.80 0.84 0.88 0.93 0.95 0.98 1.00	0.05915250	−0.07
	y:	0.00 0.02 0.04 0.05 0.06 0.07 0.09 0.10 0.12 0.14 0.15 0.16 0.18 0.20 0.21 0.22 0.24 0.26 0.27 0.28 0.30 0.32 0.33 0.34 0.36 0.38 0.39 0.40 0.42 0.43 0.45 0.46 0.48 0.49 0.51 0.53 0.54 0.55 0.57 0.58 0.60 0.61 0.63 0.65 0.66 0.68 0.69 0.71 0.72 0.73 0.75 0.77 0.78 0.80 0.84 0.88 0.93 0.95 0.98 1.00
801	x:	$x=\sqrt{z}$ , $0\le z\le 1$ , $z=0.00125i$ , $i=0,1,\cdots ,800$	0.05919402	0.004
	y:	$y=\sqrt{z}$ , $0\le z\le 1$ , $z=0.00125i$ , $i=0,1,\cdots ,800$

Figure 5. Graph of $z=\sin \left(20x\right)+\sin \left(20y\right)$ over $\left[0,1\right]\times \left[0,1\right]$ . As in the previous example, a finer discretization is required for convergence.

5. Results and Discussion

The proposed matrix-based formulation provides a compact representation of the composite Simpson’s 1/3 rule for double numerical integration over rectangular domains with prescribed non-uniformly spaced nodes. By applying Simpson’s rule over consecutive pairs of subintervals in each coordinate direction, the contributions of all subrectangles are accumulated into a fixed 3 × 3 matrix $M$ , leading to the compact expression

$I\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)$

The main structural advantage of the proposed formulation is that the resulting matrix has fixed size 3 × 3, independently of the number of discretization nodes. This provides a compact representation of the numerical integration and can reduce the auxiliary memory required for the matrix formulation. In particular, the nine entries of $M$ can be accumulated successively without requiring an additional matrix whose dimensions increase with the number of nodes. This memory reduction is relative to an implementation of the conventional composite Simpson formulation that explicitly constructs and stores the full function-value matrix. A streaming implementation of the conventional composite Simpson rule can also process function values sequentially without storing the complete function-value matrix; therefore, the proposed formulation does not claim a general memory advantage over such implementations.

However, the fixed size of the resulting matrix does not imply a reduction in the asymptotic computational complexity. Since each entry of $M$ is obtained through double summations over the subrectangles, the overall computational cost remains of order $O\left(nm\right)$ , as in the conventional composite Simpson formulation. Therefore, the principal computational advantage of the proposed approach lies in its compact representation and reduced auxiliary memory requirements rather than in a demonstrated reduction of computational time.

An important feature of the method is its ability to operate on prescribed non-uniform rectangular meshes. This is particularly relevant when function values are available at previously determined non-uniform locations. Unlike Gaussian quadrature, whose nodes are determined by the selected quadrature rule, the proposed method does not require redistribution of the prescribed nodes. Nevertheless, when only discrete experimental data are available, the midpoint values required by Simpson’s rule must also be available or obtained by an appropriate procedure.

The proposed method is complementary to Gaussian quadrature and adaptive cubature methods. Gaussian quadrature can achieve high accuracy for sufficiently smooth functions when the evaluation points can be freely selected, while adaptive cubature methods are advantageous when automatic refinement and error control are required. The present method instead emphasizes simplicity, a fixed matrix structure, and the direct treatment of prescribed non-uniform meshes.

The error analysis developed in Section 3.2 establishes fourth-order convergence for the uniform-spacing case. For general non-uniform meshes, the local mesh widths vary from one subrectangle to another, and a complete theoretical characterization of the convergence order has not been established in the present work. The fourth-order behavior observed in the numerical experiments for non-uniform meshes is therefore reported as an empirical observation and should not be interpreted as a theoretically established convergence result. The numerical results presented in this work nevertheless illustrate the applicability of the formulation to non-uniform meshes.

Overall, the proposed method provides a simple and compact matrix-based reformulation of composite Simpson’s 1/3 rule for double integration. Its main characteristics are the fixed 3 × 3 matrix representation, the ability to accommodate prescribed non-uniform nodes, and the possibility of reducing auxiliary memory requirements. These features make the method a useful complementary approach for numerical integration on non-uniform rectangular meshes.

6. Conclusions

A compact matrix-based formulation of the composite Simpson’s 1/3 rule has been developed for the numerical evaluation of double integrals over rectangular domains with prescribed non-uniformly spaced nodes. The proposed formulation reorganizes the contributions from the individual subrectangles into a fixed 3 × 3 matrix, whose entries are obtained from weighted double summations. Consequently, the numerical approximation can be expressed in the compact form

$I\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)$

The fixed-size matrix representation provides a simple and compact organization of the integration procedure, independently of the number of discretization nodes. In particular, the nine entries of the matrix can be accumulated successively, which can reduce the auxiliary memory required by the matrix representation when dealing with large datasets.

The proposed formulation preserves the fourth-order behavior of the composite Simpson’s 1/3 rule under the uniform-spacing assumptions considered in the error analysis. For prescribed non-uniform meshes, the numerical examples demonstrate that the method can provide accurate approximations while avoiding the redistribution of the available nodes. A complete theoretical analysis of the convergence order for general non-uniform meshes remains an open issue for further study.

The proposed method is complementary to Gaussian quadrature and adaptive cubature methods. While Gaussian quadrature is particularly effective when the evaluation nodes can be selected according to the quadrature rule, and adaptive cubature methods provide automatic refinement and error control, the present formulation is intended for situations in which data are available on a prescribed non-uniform rectangular mesh.

Although the proposed formulation maintains the $O\left(nm\right)$ asymptotic computational complexity associated with processing the data over the mesh, its fixed 3 × 3 structure provides a compact representation and offers the possibility of reducing auxiliary memory requirements. Therefore, the method constitutes a simple and flexible alternative for composite double numerical integration when prescribed non-uniform data are available.

Acknowledgments

The authors gratefully acknowledge the support provided by the Universidad Nacional del Santa (UNS), Chimbote, Perú, and the Universidade Federal do Pampa (Unipampa), whose institutional support contributed to the development and preparation of this work. The authors also thank the Brazilian National Council for Scientific and Technological Development (CNPq) for its financial support through the Research Productivity Fellowship awarded to Luis E. G. Armas.

Author Contributions

All authors contributed to the conception and design of the study and to the final formulation of the research. Material preparation, formula derivation, data collection, and manual numerical calculations were carried out by Lizandro B. R. Zegarra, Luis E. G. Armas, Aída B. D. Reyna, César C. Campos, Pedro G. R. Carrera, Herón J. M. Marchena, Carlos A. M. Alvarado, and María E. B. Gemin. In particular, to facilitate the extensive numerical calculations, a computer program was developed in MATLAB R2020a by Herón J. M. Marchena and Lizandro B. R. Zegarra. The first draft and the final version of the manuscript were prepared by Lizandro Baldomero Reyna Zegarra, and all authors provided comments on previous versions of the manuscript. All authors read and approved the final manuscript.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Otto, S.R. and Denier, J.P. (2005). An Introduction to Programming and Numerical Methods in MATLAB. Springer.
[2]	Süli, E. and Mayers, D.F. (2003) An Introduction to Numerical Analysis. Cambridge University Press.[CrossRef]
[3]	Nieves, A. and Dominguez, F. (2014) Metodos numericos aplicados a la ingeniería. Grupo Editorial Patria.
[4]	Grant, C. and Talvila, E. (2019) Elementary Numerical Methods for Double Integrals. arXiv: 1905.05805v1.
[5]	Lee, W.R., Wang, S. and Teo, K.L. (2000) An Optimization Approach to Numerical Integration in Two Dimensions. Applied Mathematics and Computation, 109, 205-223.[CrossRef]
[6]	Keshavarzzadeh, V., Kirby, R.M. and Narayan, A. (2018) Numerical Integration in Multiple Dimensions with Designed Quadrature. SIAM Journal on Scientific Computing, 40, A2033-A2061.[CrossRef]
[7]	Dehda, B. (2019) A New Approach to Numerical Integration Based on Coifman Wavelets. Journal of Applied Mathematics and Computational Mechanics, 18, 31-44.[CrossRef]
[8]	Asgedom, A.A. and Kefela, Y.Y. (2026) An Adaptive Hybrid Quadrature Scheme: Combining Simpson’s Rule and Gaussian Quadrature for Enhanced Numerical Integration. PLOS One, 21, e0335582.[CrossRef]
[9]	Mohanty, S. and Alahmadi, J (2026) A Generalised Cubature Rule Using Simpson Cubature Rule. Annals of Pure and Applied Mathematics, 33, 53-62.
[10]	Zhong H. and Feng X. (2024) An Efficient Implementation Algorithm for Quasimonte Carlo Approximations of High-Dimensional Integrals. arXiv: 2404.08867v1.
[11]	Yang, W.Y., Cao, W., Chung, T. and Morris, J. (2005) Applied Numerical Methods Using Matlab®. Wiley.[CrossRef]