A smoothing spline algorithm to interpolate and predict the eigenvalues of matrices extracted from the sequence of preconditioned banded symmetric Toeplitz matrices

1.   Introduction

In mathematics, linear algebra is essential for resolving real-world issues. Any numerical approach used to approximated differential equations leads to different matrices, and among them are structured, multilevel, banded matrices, and matrix-sequences play an important role. At the beginning of the 20th century, Otto Toeplitz (1881–1940), a German mathematician, introduced the notion of the Toeplitz matrix and the Toeplitz operator. After the Fourier transformation was discovered, in order to determine the spectral theory of the Toeplitz matrix associated with the Laurent operators $L = (\Phi_{j-i})_{i, j\in Z}$ on the Lebesgue space $ł^2(z)$, Otto Toeplitz employed the multiplication operation. This resulted in the following:

Stated otherwise, this defines the Toeplitz matrix of size $n$ with a fixed entry along each diagonal as

Grenander and Szego ^[4] state that the $n^{\text{th}}$ Toeplitz matrix generated by $\Phi$, for a complex valued Lebesgue integrable function $\Phi:[-\pi, \pi]\to\mathbb C$, is given by

where the Fourier coefficients of $\Phi$ are denoted by the quantities $\hat \Phi_k$, i.e.,

The Toeplitz sequence created by $\Phi$, also known as the generating function of $\{T_n(\Phi)\}_n$, is denoted as $\{T_n(\Phi)\}_n$. The spectral properties of all the matrices $T_n(\Phi)$ where $\Phi$ is real-valued are well-known, ranging from the eigenvalue localization to the Weyl sense asymptotic spectral distribution. Specifically, if $\Phi$ is real-valued, $\Phi$ represents $\{T_n(\Phi) \}$ spectrally; see, for example, ^[18] and the references therein for further details.

More specifically, any eigenvalue of $T_n(\Phi)$ with $\Phi$ being real-valued and not identically constant falls in $(m_\Phi, M_\Phi)$, where $m_\Phi$ and $M_\Phi$ represent the infimum and supremum of $\Phi$, respectively. In the event that $\Phi$ is constant, there is no problem; if $\Phi = m$ then $T_n(\Phi) = m I_n$, where $I_n$ represents the identity of size $n$. Therefore, $T_n(\Phi)$ is Hermitian positive definite if $M_\Phi > 0$ and $\Phi$ is non-negative almost everywhere. The eigenvalues of banded symmetric Toeplitz (BST) matrices asymptotically expand when the generating symbol is a real-valued cosine trigonometry polynomial (RCTP) in the form specified by Ekstrm et al. ^[7]

Over the interval $[0, \pi]$, the RCTP is either monotonically increasing or decreasing. The $n\times n$ real banded symmetric Toeplitz (BST) matrix is generated by $\Omega$. We focus on the following assumptions in this paper: $T_n(\Omega)$, $T_n(\Phi)$ are both real symmetric Toeplitz matrices, given

● $T_n(\Phi)$ is a positive definite matrix.

● We take ${\mathcal P}_n(\Omega, \Phi) = T_n^{-1}(\Omega) T_n(\Phi)$ as the 'preconditioned' matrix and define $r = \Phi/\Omega$.

The Toeplitz matrix of order $n \times n$ with bandwidth $2m+1$ was generated by the symbol $p\in \{\Omega, \Phi\}$, where $m\in \{m_1, m_2\}$ (notice that $m = m_1$ when $p = \Omega$ and $m = m_2$ when $p = \Phi$). Consequently, the matrix is

The quick solution of large Toeplitz linear systems requires matrices of the form ${\mathcal P}_n(\Omega, \Phi)$ (in relation to the preconditioned conjugate gradient method). Their spectral characteristics have been thoroughly investigated. To be more precise, if $r = m$, then ${\mathcal P}_n(\Omega, \Phi) = r I_n$; on the other hand, all eigenvalues of ${\mathcal P}_n(\Omega, \Phi)$ belong to $(m_r, M_r)$ if $m_r < M_r$ (see, for example, ^[21]) and $r = \Phi/\Omega$ is the Weyl sense spectral distribution for the complete sequence $\{{\mathcal P}_n(\Omega, \Phi)\}_n$ ^[24]. If a function increases or decreases over the interval $[0, \pi]$, it is considered monotone in our context. Assuming that $r = \Phi/\Omega$ is monotone, we experimentally demonstrate in this article that the following asymptotic expansion holds ^[5] for every $n$, every integer $\alpha\ge0$, and every $k = 1, \ldots, n$.

where

● eigenvalues of ${\mathcal P}_n(\Omega, \Phi)$ are arranged in either non-increasing or non-decreasing order, contingent on whether $r$ is increasing or decreasing;

● a series of functions from $[0, \pi]$ to $\mathbb R$ that depends entirely on $r$ is denoted as $\{ \zeta_j\}_{j = 1, 2, \ldots}$;

● $\Lambda = \frac{1}{n+1}$ and $\Theta_{k, n} = \frac{k\pi}{n+1} = k\pi \Lambda$;

● the remainder (error) is $\epsilon_{j, n, \alpha} = O({\Lambda}^{\alpha+1})$. This results in an inequality for some constant $C_\alpha$ that depends only on $\alpha$ and $r$: $|\epsilon_{k, n, \alpha}|\le C_\alpha {\Lambda}^{\alpha+1}$.

If the RCTP $\Phi$ is monotone and satisfies other conditions, like $\Phi'(\Theta)\ne0$ and $\Phi''(\Theta)\ne0$ for $\Theta\in\{0, \pi\}$, the result is proven in ^[14,16,17]. This is the pure Toeplitz case, meaning that ${\mathcal P}_n(\Omega, \Phi) = T_n(\Phi)$ and $r = \Phi$ are equivalent for $\Omega = 1$. Particularly interesting are

that appear when discretizing differential equations. Sadly, if $q\ge2$, the condition that $\Phi''(0)\ne0$ is not satisfied. Even in this "degenerate case" the higher order approximation (1.3) holds, according to numerical evidence in ^[22]. Here, given the preconditioned matrices ${\mathcal P}_n(\Omega, \Phi)$, we numerically prove the same. Theoretically, the demonstration of the aforementioned conjecture when $\alpha = 0$ complements the numerical testing.

The authors of ^[22] utilized the asymptotic expansion (1.3) to approximate $\upsilon_k(T_n(\Omega))$ for very large $n$, given that the values for fairly sized $n_1, \ldots, n_s$ with $\Theta_{k_1, n_1} = \cdots = \Theta_{k_s, n_s} = \Theta_{k, n}$, $s\ge 2$ are available. This paper's second goal is to implement this theory and provide evidence for it through numerical experiments and a suitable spline error analysis. Specifically, we provide a method that can compute $\upsilon_k({\mathcal P}_n(\Omega, \Phi)$ with low computational cost and high degree of accuracy. The method is comparable to the extrapolation process used in Romberg integration to produce high-precision integral estimates from a small number of imprecise trapezoidal approximations. This is where the Euler-Maclaurin formula and the asymptotic expansion (1.3) come into play.

2.   Algorithm

Throughout the entire work, we make use of the following parameters: Step size $\Lambda = 1/(n+1)$, for certain positive integers $n\in \mathbf{N}$, and the grid points $\Theta = k\pi \Lambda, \, k = 1, 2, 3\cdots, n$. We make some assumptions and employ expansion (1.3) throughout this section:

● The RCTPs are $\Omega, \, \Phi, \, r$, and they increase monotonically.

● Suppose that $\alpha$ and $\, n(1), n \in \mathbf{N}$, are fixed parameters.

● Either $n_j = m^{j-1}(n(1)+1)-1$ or $n_j$ are selected at random.

● Give definitions to the PBST matrix sequences ${\mathcal P}_{n_j}(\Omega, \Phi)$.

● Determine the eigenvalues of each and every PBST matrix, ${\mathcal P}_{n_j}(\Omega, \Phi)$.

This algorithm is intended to compute the eigenvalues of the ${\mathcal P}_n(\Omega, \Phi)$ when $n$ is very large, which are very difficult to compute by the standard eigensolver (e.g., MATLAB function, 'eig'). Therefore, in order to expedite computation, we compute the eigenvalues of the ${\mathcal P}_{n_1}, \, {\mathcal P}_{n_2}, \cdots, {\mathcal P}_{n_s}$, for small $n_j$ by the eigensolver. Using these eigenvalues of the small size PBST matrices, we can evaluate the eigenvalues of the large size PBST matrix.

Interpolation and extrapolation

When $\Theta_{j, n_1} = \Theta_{j, n_s}$ is fixed, we can compute the eigenvalues of all ${\mathcal P}_{n_j}(\Omega, \Phi)$ and set them in Table 1 against ${\Lambda}_j$ as

Table 1.  Eigenvalues of ${\mathcal P}_{n_j}$ matrices according to $s$.

${\Lambda}_1$	$\upsilon_{1, 1}$	$\upsilon_{2, 1}$	$\cdots$	$\upsilon_{n_1, 1}$
${\Lambda}_2$	$\upsilon_{1, 2}$	$\upsilon_{2, 2}$	$\cdots$	$\upsilon_{n_2, 2}$
${\Lambda}_3$	$\upsilon_{1, 3}$	$\upsilon_{2, 3}$	$\cdots$	$\upsilon_{n_3, 3}$
$\vdots$	$\vdots$	$\vdots$	$\cdots$	$\vdots$
${\Lambda}_{s}$	$\upsilon_{1, s}$	$\upsilon_{2, s}$	$\cdots$	$\upsilon_{n_s, s}$

---

 | Show Table

DownLoad: CSV

Every collection of eigenvalues needs to have a grid attached to it. In order to define the grid, we define $\Theta_{n_j, k} = k\, {\Lambda}_j\, \pi$ for $k = 1, 2, 3, \cdots, n_j$, and ${\Lambda}_j = \frac{1}{n_j+1}$. For $k = 1, 2, 3, \cdots n_\beta$, the $\Theta$-grid of the ${\mathcal P}_{n_\beta}(\Omega, \Phi)$ is $\Theta_{n_\beta, j} = j\, {\Lambda}\, \pi$. For $j = 1, 2, 3, \cdots, s$, the number of eigenvalues in ${\mathcal P}_{n_\beta}(\Omega, \Phi)$ is significantly more than that in ${\mathcal P}_{n_j}(\Omega, \Phi)$. For each set of eigenvalues of size $n_j$, $j = 1, 2, 3, \cdots, s$ attached to $\Theta$-grid $\Theta_{n_j, k}$, $k = 1, 2, 3, \cdots, n_\beta$, we interpolate and extrapolate the $n_\beta$-eigenvalues before doing the extrapolation $k = n_j, \cdots, 1, 2, 3$. We apply a higher-order spline-curve fitting for all $\Theta$-grids in order to attain the higher order of precision in the interpolation and extrapolation of eigenvalues. For the interior $\Theta$-grid node, the suggested higher-order spline fitting performs well; nevertheless, performance degrades toward the grid's edge. We apply lower-order shape-preserving spline curve fitting to solve this issue.

After all $\Theta$-grids have been interpolated and extrapolated, we perform one last extrapolation to obtain the eigenvalues of ${\mathcal P}_{n_\beta}(\Omega, \Phi)$. We have all of the interpolated and extrapolated eigenvalues of the PBST matrices ${\mathcal P}_{n_j}(\Omega, \Phi)$. We can now compute the $\tilde {\zeta}_k(\Theta_{n_j, j})$ using these eigenvalues. We extend the expansion (1.3) for $\alpha = 4$ in order to achieve this.

Set $n_j$, $j\in\{1, 2, 3, 4 \}$ for any positive integers $m$ and $n_j$, then satisfy $n_j = m^{j-1}(n(1)+1)-1$ or pick $n_j$ randomly. The expansion (2.1) for the four matrices of $n_j$ is as follows:

Observe that $\Theta_{k_i, n_i} = \Theta_{k_1, n_1} = \bar{\Theta}$ for $k_1\in\{1, 2, \cdots, n_1\}$ where ${\Lambda}_i = \frac{1}{n_i+1}$. The numerical estimate of $\zeta_j(\bar{\Theta})$ for $j\in\{1, 2, 3, 4\}$ is of interest to us. The set of Eq (2.2) becomes

For $k_1 \in 1, 2, \cdots, n_1$, we solve the aforementioned system of linear equations to obtain $\tilde {\zeta}_k(\bar {\Theta})$. The big size $n_{\beta}$ matrix's eigenvalues are estimated using the computed $\tilde {\zeta}_j$ by utilizing

References for the error bounds of \(\zeta_k\) in the asymptotic expansion are provided in ^[5].

3.   Description of the algorithm

Finding all of the eigenvalues of the large-size PBST matrices is the primary goal of the numerical algorithm. Using conventional eigenvalue solvers such as the 'eig' function in MATLAB, we can obtain the spectrum information of smaller PBST matrices. The generating symbol is the quotient of RCTPs since we are working with PBST matrices. In order for the suggested technique to achieve high accuracy in the numerically computed eigenvalues, the generating symbol's monotonicity is essential. The algorithm remains functional even when dealing with non-monotonically generated symbols; however, the non-monotonic part's numerical accuracy suffers as a result.

Let $\{{\mathcal P}_n(\Omega, \Phi) = T_n^{-1}(\Phi) T_n(\Omega)\}$ be a sequence of preconditioned symmetric banded Toeplitz matrices. For $k = n_1, n_2, n_3, \cdots n_s$, we select a collection of ${\mathcal P}_k(\Omega, \Phi)$. The total number of small size ${\mathcal P}_k(\Omega, \Phi)$ is denoted by $s$ in this case. Furthermore, the $n_j$ are ordered as follows: $n_1 < n_2 < n_3 < \cdots < n_s$. We do not impose any stringent requirements on selecting $n_j$'s, but careful selection can raise the accuracy of our numerical findings. We might calculate the eigenvalues of ${\mathcal P}_{n_j}(\Omega, \Phi)$ for each $n_j$. For $i = 1, 2, 3, \cdots, s$, we have $s$-sets of eigenvalues with sizes $n_j$. Assume that $n_\beta$ is a huge number and that we are interested in calculating the spectrum of a large size ${\mathcal P}_{n_\beta}(\Omega, \Phi)$. Every collection of eigenvalues needs to have a grid attached to it. In order to define the grid, we define $\Theta_{n_j, k} = k\, {\Lambda}_j\, \pi$ for $k = 1, 2, 3, \cdots, n_j$, and ${\Lambda}_j = \frac{1}{n_j+1}$. For $k = 1, 2, 3, \cdots n_\beta$, the $\Theta$-grid of the ${\mathcal P}_{n_\beta}(\Omega, \Phi)$ is $\Theta_{n_\beta, j} = j\, {\Lambda}\, \pi$.

For $j = 1, 2, 3, \cdots, s$, the number of eigenvalues in ${\mathcal P}_{n_\beta}(\Omega, \Phi)$ is significantly more than that in ${\mathcal P}_{n_j}(\Omega, \Phi)$. For each set of eigenvalues of size $n_j$, $j = 1, 2, 3, \cdots, s$ attached to $\Theta$-grid $\Theta_{n_j, k}$, $k = 1, 2, 3, \cdots, n_\beta$, we interpolate and extrapolate the $n_\beta$-eigenvalues before doing the extrapolation. We employ higher-order spline-curve fitting across all \(\Theta\)-grids to achieve greater precision in interpolating and extrapolating eigenvalues. For nodes within the \(\Theta\)-grid, this advanced spline fitting shows commendable performance; however, its efficacy diminishes near the grid's boundaries. To address this, we resort to lower-order, shape-preserving spline fitting. Once interpolation and extrapolation are completed for all \(\Theta\)-grids, a final extrapolation step is undertaken to ascertain the eigenvalues of \({\mathcal P}_{n_\beta}(\Omega, \Phi)\).

It is worth mentioning that Bogoya et al. ^[6] recently published an article introducing a new algorithm that can accurately interpolate and extrapolate the spectrum of preconditioned Toeplitz matrix sequences at machine level accuracy. The approach described by Bogoya et al. ^[6] differs from the current idea, primarily in that it requires the generating symbol to be directly inverted. Alternatively, one can solve a nonlinear equation to get the inversely projected values; however, our current methodology does not require this step.

An algorithm for extracting particular eigenvalues from huge preconditioned matrices generated from a sequence of preconditioned Toeplitz matrices was published by Fayyaz et al. ^[5]. Our current work is novel since we can extrapolate and interpolate the whole spectrum. One limitation of our strategy is that it necessitates the use of two different kinds of smoothing-spline methods, particularly in proximity to boundaries where a shape-preserving smoothing-spline approach is required.

4.   Numerical testing

Here we perform numerical testing on a collection of problems taken from the paper of Ekström et al. ^[7]. We assume $\Theta \in [0, \pi]$ and $\Omega(\Theta) > 0, \, \forall \, \Theta \in [0, \pi]$, $\Phi''(0)\neq 0$, for all of our numerical testing.

Example 1. Consider the monotonic rising functions $\Omega$, $\Phi$, and $r$, which are defined as

For $n = 5000, \, n(1) = 10$, and $s = 7$, we wish to approximate the eigenvalues of the $\upsilon_n = T_n(\Omega)^{-1}T_n{\Phi} = T_n(r)$. The error between eigenvalues $\upsilon_n$ and an approximated eigenvalues $\tilde{\upsilon_n}$ is displayed in Figure 1, when compared to $\Theta_{j, n}, \, j = 1, 2, 3, \cdots, n$. The graph illustrates the compression between two algorithms, the first of which presented in ^[7]. Compared to the Erik Algorithm ^[7], our proposed technique is faster and more accurate, as shown by the Figure 1.

Figure 1.  $log10$(Error), in the case $n = 5000$, $n(1) = 10$, and $s = 7$ for Example 1.

DownLoad: Full-Size Img PowerPoint

Example 2. Let the monotonic rising functions $\Omega$, $\Phi$, and $r$ be defined as

For $n = 10000, \, n(1) = 10$, and $s = 7$, we wish to approximate the eigenvalues of the $\upsilon_n = T_n(\Omega)^{-1}T_n{\Phi} = T_n(r)$. The error between the approximated eigenvalues $\tilde{\upsilon_n}$ and the eigenvalues $\upsilon_n$ is displayed in Figure 2 in relation to $\Theta_{j, n}, \, j = 1, 2, 3, \cdots, n$. The graph indicates that our algorithm's accuracy increases as the target matrix's size increases, with values above machine zero.

Figure 2.  $log10$(Error), in the case $n = 10000$, $n(1) = 10$, and $s = 7$ for Example 2.

DownLoad: Full-Size Img PowerPoint

Example 3. Let the monotonic rising functions $\Omega$, $\Phi$, and $r$ be described as

For $n = 10000, \, n(1) = 10$, and $s = 5$, we wish to approximate the eigenvalues of the $\upsilon_n = T_n(\Omega)^{-1}T_n{\Phi} = T_n(r)$. The error between the eigenvalues $\upsilon_n$ and the approximated eigenvalues $\tilde{\upsilon_n}$ versus $\Theta_{j, n}, \, j = 1, 2, 3, \cdots, n$, is shown in Figure 3. When we increase the value of $s$, it means the sequence of the PBST matrices increases. We can observe this immediately.

Figure 3.  $log10$(Error), in the case $n = 10000$, $n(1) = 10$, and $s = 5$ for Example 3.

DownLoad: Full-Size Img PowerPoint

Example 4. Let the monotonic rising functions $\Omega$, $\Phi$, and $r$ be described as

Next, the eigenvalues of $\upsilon_n = T_n(\Omega)^{-1}T_n{\Phi} = T_n(r)$ are approximated for $n = 10000, \, n(1) = 10$, and $s = 7$. The error between the eigenvalues $\upsilon_n$ and the approximated eigenvalues $\tilde{\upsilon_n}$ versus $\Theta_{j, n}, \, j = 1, 2, 3, \cdots, n$, is virtually minimum from Example $4$ in ^[7], as we can see in Figure 4. Our algorithm errors are quickly evaluated and range from $10^{-9}$ to over $10^{-16}$.

Figure 4.  $log10$(Error), in the case $n = 10000$, $n(1) = 10$, and $s = 7$ for Example 4.

DownLoad: Full-Size Img PowerPoint

Example 5. Consider the monotonic rising functions $\Omega$, $\Phi$, and $r$, which are defined as

We wish to approximate the eigenvalues of the $\upsilon_n = T_n(\Omega)^{-1}T_n{\Phi} = T_n(r)$ for $n = 10000, \, n(1) = 25$, and $s = 7$. The error between the estimated eigenvalues $\tilde{\upsilon_n}$ and the eigenvalues $\upsilon_n$ is displayed in Figure 5 in relation to $\Theta_{j, n}, \, j = 1, 2, 3, \cdots, n$. The lowest error is obtained for $\Theta_{j, n}$ when $j \in [500, 4500]$.

Figure 5.  $log10$(Error), in the case $n = 10000$, $n(1) = 25$, and $s = 7$ for Example 5.

DownLoad: Full-Size Img PowerPoint

Example 6. Let the monotonic rising functions $\Omega$, $\Phi$, and $r$ be defined as

We wish to approximate the eigenvalues of the $\upsilon_n = T_n(\Omega)^{-1}T_n{\Phi} = T_n(r)$ where $r = \Phi/\Omega$ for $n = 5000, \, n(1) = 50$, and $s = 4$. The error between the estimated eigenvalues $\tilde{\upsilon_n}$ and the eigenvalues $\upsilon_n$ is displayed in Figure 6 in relation to $\Theta_{j, n}, \, j = 1, 2, 3, \cdots, n$. The lowest error is obtained for $\Theta_{j, n}$ when $j \in [500, 4500]$. The discrepancy between the eigenvalues can be observed by varying the parameters and producing functions that are more adequate for all instances, as demonstrated in Examples $1$ to $5$ when $\Omega(\Theta)$ is similar. These codes are found in ^[7]. As we can see in Figure 6, there are more errors than machine zeros in case 6, where $\Omega(\Theta)$ is not identical. When both $\Phi(\Theta)$ and $\Omega(\Theta)$ are not equal, our algorithm performed most accurately.

Figure 6.  $log10$(Error), in the case $n = 5000$, $n(1) = 50$, and $s = 4$ for Example 6.

DownLoad: Full-Size Img PowerPoint

Table 2.  Example 6. CPU time, in the case $\Omega(\Theta) = 3+2 \cos(\Theta), \, \Phi(\Theta) = 2- \cos(\Theta)- \cos(2\Theta)$.

Parameters	Erik Algorithm [7](Time)	Proposed Algorithm(Time)
$s=4, n=5000, n(1)= 50$	$13.3093$	$0.8154$
$s=4, n=5000, n(1)= 100$	$12.8770$	$1.6483$
$s=4, n=5000, n(1)= 200$	$18.8837$	$3.9735$
$s=4, n=5000, n(1)= 400$	$25.7752$	$14.0923$
MATLAB's eig function	$44.9260$	$41.1051$

---

 | Show Table

DownLoad: CSV

Example 7. Let the monotonic rising functions $\Omega$, $\Phi$, and $r$ be defined as

We wish to approximate the eigenvalues of the $\upsilon_n = T_n(\Omega)^{-1}T_n{\Phi} = T_n(r)$ for $n = 5000, \, n(1) = 50$, and $s = 4$. In Figure 7, the error between the eigenvalues $\upsilon_n$ and an approximated eigenvalues $\tilde{\upsilon_n}$ versus $\Theta_{n_{\beta}, j}, \, j = 1, 2, 3, \cdots, n_{\beta}$ shows the minimum error inside of the boundary, and the maximum near the boundary, but we achieve the minimum error, more specifically the machine zero, demonstrating the high accuracy and speed of our proposed algorithm in comparison to that of Ekström ^[7].

Figure 7.  $log10$(Error), in the case $n = 5000$, $n(1) = 50$, and $s = 4$ for Example 7.

DownLoad: Full-Size Img PowerPoint

5.   Conclusions

In this work, we present an algorithm to calculate the eigenvalues of matrices extracted from a sequence of large-scale Toeplitz matrices that are preconditioned banded symmetric (PBST). When compared to existing eigenvalue solvers, computing the eigenvalues of large-size PBST matrices should be more economical. We collect all the spectrum information of small-size Toeplitz PBST matrices while taking the algorithm's computing cost into consideration. The reciprocal sizes of smaller PBST matrices can be used to sort their spectral information. Depending on the size of the PBST matrix, the monotonic spectrum arrangement can be arranged in bijection with an appropriate grid. For the large-size grid connected to the large-size PBST, we are able to interpolate and extrapolate the eigenvalue of the small-size PBST. Now, all of the small PBSTs have the same size spectrum. Finally, we use the extrapolation method---in our case, higher-order spline fitting---to compute the spectrum of the large-size PBST.

Author contributions

Salima Kouser: Methodology, Writing- Original draft preparation. Shafiq Ur Rehman: Conceptualization, Supervision. Mabkhoot Alsaiari: Investigation, Formal analysis. Fayyaz Ahmad: Software, Validation. Mohammed Jalalah: Visualization, Resources, Funding acquisition. Farid A. Harraz: Data curation, Writing - Review & editing. Muhammad Akram: Project administration, Writing - Review & editing.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgment

The authors are grateful to the Deanship of Scientific Research at Najran University for funding this work under the Research Groups Funding Programme grant code (NU/RG/SERC/12/23).

Conflict of interest

The authors declare that there is no conflict of interest regarding the publication of this paper.