Symmetry, Hopf bifurcation, and offset boosting in a novel chameleon system

Jie Liu; Bo Sang; Lihua Fan; Chun Wang; Xueqing Liu; Ning Wang; Irfan Ahmad; Jie Liu; Bo Sang; Lihua Fan; Chun Wang; Xueqing Liu; Ning Wang; Irfan Ahmad

doi:10.3934/math.2025225

AIMS Mathematics

2025, Volume 10, Issue 3: 4915-4937. doi: 10.3934/math.2025225

Previous Article Next Article

Research article Special Issues

Symmetry, Hopf bifurcation, and offset boosting in a novel chameleon system

1.
School of Mathematical Sciences, Liaocheng University, Liaocheng 252059, Shandong, China
2.
School of Microelectronics and Control Engineering, Changzhou University, Changzhou 213159, Jiangsu, China
3.
Faculty of Science and Engineering, Maynooth International Engineering College, Maynooth University, Maynooth, Co. Kildare W23, Ireland

Received: 23 December 2024 Revised: 02 February 2025 Accepted: 17 February 2025 Published: 06 March 2025
MSC : 37D45, 37G35

Chameleon systems are dynamical systems that exhibit either self-excited or hidden oscillations depending on the parameter values. This paper presents a comprehensive investigation of a quadratic chameleon system, including an analysis of its symmetry, dissipation, local stability, Hopf bifurcation, and various chaotic dynamics as the control parameters $(\mu, a, c)$ vary. Here, $\mu$ serves as the dissipation parameter in the $y$ direction. Bifurcation analysis for four scenarios with $\mu = 0$ was performed, revealing the emergence of various dynamical phenomena under different parameter settings. Offset boosting means introducing a constant into one of the state variables of the system for boosting the variable to a different level. Additionally, hidden chaotic bistability with offset boosting was exhibited by varying $\mu$ . The parameter $\mu$ serves as both the Hopf bifurcation parameter and the offset boosting parameter, while the other parameters $(a, c)$ also play critical roles as control parameters, resulting in period-doubling routes to self-excited or hidden chaotic attractors. These findings enrich our understanding of nonlinear dynamics in quadratic chameleon systems.

Keywords:

Citation: Jie Liu, Bo Sang, Lihua Fan, Chun Wang, Xueqing Liu, Ning Wang, Irfan Ahmad. Symmetry, Hopf bifurcation, and offset boosting in a novel chameleon system[J]. AIMS Mathematics, 2025, 10(3): 4915-4937. doi: 10.3934/math.2025225

Related Papers:

[1]	Ali H. Tedjani, Sharifah E. Alhazmi, Samer S. Ezz-Eldien . An operational approach for one- and two-dimension high-order multi-pantograph Volterra integro-differential equation. AIMS Mathematics, 2025, 10(4): 9274-9294. doi: 10.3934/math.2025426
[2]	Ahmed M. Rajab, Saeed Pishbin, Javad Shokri . Analyzing the structure of solutions for weakly singular integro-differential equations with partial derivatives. AIMS Mathematics, 2024, 9(9): 23182-23196. doi: 10.3934/math.20241127
[3]	Obaid Algahtani, M. A. Abdelkawy, António M. Lopes . A pseudo-spectral scheme for variable order fractional stochastic Volterra integro-differential equations. AIMS Mathematics, 2022, 7(8): 15453-15470. doi: 10.3934/math.2022846
[4]	A. H. Tedjani, A. Z. Amin, Abdel-Haleem Abdel-Aty, M. A. Abdelkawy, Mona Mahmoud . Legendre spectral collocation method for solving nonlinear fractional Fredholm integro-differential equations with convergence analysis. AIMS Mathematics, 2024, 9(4): 7973-8000. doi: 10.3934/math.2024388
[5]	Chuanli Wang, Biyun Chen . An $hp$ -version spectral collocation method for fractional Volterra integro-differential equations with weakly singular kernels. AIMS Mathematics, 2023, 8(8): 19816-19841. doi: 10.3934/math.20231010
[6]	Mohammed A. Almalahi, Satish K. Panchal, Fahd Jarad, Mohammed S. Abdo, Kamal Shah, Thabet Abdeljawad . Qualitative analysis of a fuzzy Volterra-Fredholm integrodifferential equation with an Atangana-Baleanu fractional derivative. AIMS Mathematics, 2022, 7(9): 15994-16016. doi: 10.3934/math.2022876
[7]	Chuanhua Wu, Ziqiang Wang . The spectral collocation method for solving a fractional integro-differential equation. AIMS Mathematics, 2022, 7(6): 9577-9587. doi: 10.3934/math.2022532
[8]	Khadijeh Sadri, David Amilo, Kamyar Hosseini, Evren Hinçal, Aly R. Seadawy . A tau-Gegenbauer spectral approach for systems of fractional integro-differential equations with the error analysis. AIMS Mathematics, 2024, 9(2): 3850-3880. doi: 10.3934/math.2024190
[9]	Mian Bahadur Zada, Muhammad Sarwar, Reny George, Zoran D. Mitrović . Darbo-Type $\mathcal{Z}_{\rm{m}}$ and $\mathcal{L}_{\rm{m}}$ contractions and its applications to Caputo fractional integro-differential equations. AIMS Mathematics, 2021, 6(6): 6340-6355. doi: 10.3934/math.2021372
[10]	Hawsar Ali Hama Rashid, Mudhafar Fattah Hama . Approximate solutions for a class of nonlinear Volterra-Fredholm integro-differential equations under Dirichlet boundary conditions. AIMS Mathematics, 2023, 8(1): 463-483. doi: 10.3934/math.2023022

Abstract

1. Introduction

The goal of this paper is to develop Legendre spectral tau approaches for the following high-order Pantograph Volterra-Fredholm integro-differential equations (P-VF-IDEs):

● One-dimensional case:

$\begin{equation} \begin{split} \sum\limits_{k = 0}^{n}b_k y^{(k)}(x) = g(x)+\sum\limits_{k = 0}^{n}c_k y^{(k)}(p_k x)+\sum\limits_{k = 0}^{n}\int_{0}^{p_k x}d_k y^{(k)}(s)ds+\int_{0}^{x}e y(s)ds+\int_{0}^{\sigma}f y(s)ds, \end{split} \end{equation}$

(1.1)

subject to

$\begin{split} \frac{d^{k} y(0) }{d x^{k}} = \varrho_k, \qquad \quad k = 0, 1, \cdots, n-1. \end{split}$

● Two-dimensional case:

$\begin{equation} \begin{split} &\sum\limits_{k_1 = 0}^{n_1}\sum\limits_{k_2 = 0}^{n_2}b_{k_1, k_2}\frac{\partial^{k_1+k_2} y(x_1, x_2) }{\partial x_1^{k_1}\partial x_2^{k_2}} = g(x_1, x_2)+\int_{0}^{\sigma_2}\int_{0}^{\sigma_1}c y(s_1, s_2) ds_1ds_2 + \int_{0}^{x_2}\int_{0}^{x_1}d y(s_1, s_2) ds_1ds_2 \\ & + \sum\limits_{k_1 = 0}^{n_1}\sum\limits_{k_2 = 0}^{n_2}\int_{0}^{p_{1} x_1}\int_{0}^{p_{2} x_2} e_{k_1, k_2} \frac{\partial^{k_1+k_2} y(s_1, s_2)}{\partial x_1^{k_1}\partial x_2^{k_2}} ds_1ds_2+ \sum\limits_{k_1 = 0}^{n_1}\sum\limits_{k_2 = 0}^{n_2}f_{k_1, k_2}\frac{\partial^{k_1+k_2} y(p_1x_1, p_2x_2)}{\partial x_1^{k_1}\partial x_2^{k_2}}, \end{split} \end{equation}$

(1.2)

subject to

$\begin{split} \frac{\partial^{k_1} y(0, x_2) }{\partial x_1^{k_1}} = \varrho_{1, k_1}(x_2), \qquad \quad k_1 = 0, 1, \cdots, n_1-1, \quad 0\leq x_2 \leq \sigma_2, \\ \frac{\partial^{k_2} y(x_1, 0) }{\partial x_2^{k_2}} = \varrho_{2, k_2}(x_1), \qquad \quad k_2 = 0, 1, \cdots, n_2-1, \quad 0\leq x_1 \leq \sigma_1.\\ \end{split}$

The problem under consideration involves several key parameters: the coefficients $b_k, c_k, d_k, e, f$ (1D case) and $b_{k_1, k_2}, e_{k_1, k_2}, f_{k_1, k_2}, c, d$ (2D case), which weight the contributions of derivatives, delay terms, and integral operators; the pantograph-type proportional delays $p_k \in (0, 1)$ (1D) and $p_1, p_2 \in (0, 1)$ (2D), introducing proportional delays such as $y(p_k x)$ ; the fixed integration limits $\sigma, \sigma_1, \sigma_2 > 0$ , defining the domains of Fredholm integrals; and the initial conditions $\varrho_k$ (1D) and $\varrho_{1, k_1}(x_2), \varrho_{2, k_2}(x_1)$ (2D), which prescribe initial values for $y$ and its derivatives. Special cases demonstrating the problem's generality include lower-order systems where $n = 0$ reduces the 1D equation to a delay Volterra-Fredholm integral equation; non-delayed systems with $p_k = 1$ , recovering classical integro-differential equations; pure Volterra dynamics achieved by omitting $f$ (1D) or the Fredholm term (2D) to isolate Volterra behavior, or removing $e$ (1D) or the Volterra term (2D) to focus on Fredholm interactions; single-integral equations formed by discarding mixed integrals; and pantograph-free cases with $c_k = 0$ and $f_{k_1, k_2} = 0$ . These examples illustrate the framework's flexibility in unifying several problems while addressing novel combinations of high-order derivatives, multi-delays, and mixed integrals.

High-order P-VF-IDEs represent a class of mathematical models that integrate effects of time delays proportional to the independent variable, known as pantograph terms, with both Volterra and Fredholm-type integro-differential equations. These equations arise in various fields, including control theory, population dynamics, and systems biology, where they model processes influenced by historical states and time-dependent interactions ^[1,2]. A pantograph is a mechanical device that collects electric current from overhead wires to supply power to electric trains and trams. Its name originates from a historical mechanical device used for copying and scaling drawings based on similar principles of motion transfer and scaling. In electric transportation systems, the pantograph ensures a consistent connection with overhead conductors, facilitating a continuous and stable flow of electricity to the vehicle. The term also originates from early drafting instruments used for geometric replication and proportional scaling ^[3,4]. Numerical approximations of P-VF-IDEs face significant challenges due to the integration of proportional delays and mixed integral terms. These complexities necessitate the use of adaptive numerical methods to achieve high accuracy. Common numerical techniques include the collocation approach ^[5,6], the direct operational tau approach ^[7], the spectral tau approach ^[8], the Taylor polynomial scheme ^[9], the Adomian decomposition method ^[10], and the Laguerre matrix scheme ^[11].

Over the years, spectral methods have proven highly effective for approximating solutions to differential and integral equations, particularly when the data and solutions are smooth. These methods use global polynomial approximations, making them especially suitable for problems requiring high accuracy. Since the solution of a delay differential equation (DDE) globally depends on its history due to the delay variable, a global spectral method presents a good approach for solving DDEs by capturing the entire solution behavior with high precision across the domain. For instance, Rahimkhani et al. ^[3] proposed a numerical method for solving fractional pantograph differential equations using the Müntz-Legendre wavelet approximation. This approach emphasized the operational matrix of fractional-order integration to reduce the differential equation to a system of algebraic equations, facilitating efficient numerical solutions. Ezz-Eldien and Doha ^[12] proposed a numerical method for solving pantograph Volterra integro-differential equations using a Chebyshev collocation technique. Yang and Hou ^[14] extended the Jacobi spectral method to nonlinear fractional pantograph differential equations by transforming them into Volterra-Fredholm integral equations. Ghomanjani et al. ^[15] introduced a Bezier curve-based approach for solving Volterra delay-integro-differential equations and linear systems with distributed input delay and saturation. Their method involved a two-step approximation strategy, where the time interval was divided into subintervals, and piecewise Bezier polynomials were applied. Huang et al. ^[16] utilized the discontinuous Galerkin method for the numerical treatment of pantograph-type differential equations. Qin et al. ^[17] developed an $hp$ -version fractional spectral collocation method for weakly singular Volterra integro-differential equations with vanishing delay, demonstrating that approximations in a fractional polynomial space can achieve exponentially decreasing errors even in the presence of weak singularities. Zaky et al. ^{[18,19,20,21]} formulated and analyzed spectral collocation schemes for solving broad classes of linear/nonlinear integral/differential equations. Rahimkhani et al. ^[22] introduced the Hahn wavelet method for fractional-order integro-differential equations, transforming them into integer-order forms using the Laplace transform. This approach allows for continuous and differentiable solutions without requiring operational matrices.

This study enhances the numerical treatment of high-order P-VF-IDEs by addressing a literature gap through the following contributions. The Legendre spectral tau method is introduced to simultaneously manage high-order derivatives, multiple pantograph delays, and mixed Volterra-Fredholm integrals within a unified framework. While existing methods, such as Chebyshev collocation ^[12,13], Jacobi spectral techniques ^[14], or Müntz-Legendre wavelets ^[3], focus on specific aspects (e.g., Volterra-only integrals or one-dimensional cases), The proposed framework unifies these elements, tackling equations previously considered intractable due to their combined complexity. By integrating high-order derivatives, multiple delays, and mixed integrals, this approach goes beyond the disjointed methods found in current literature, establishing a basis for future advancements in fractional or partial P-VF-IDEs.

This paper is organized as follows: In the next section, we provide some necessary properties of Legendre polynomial approximations. In Section 3, we construct the Legendre tau scheme for solving the one-dimensional high-order P-VF-IDEs. In Section 4, we extend the numerical scheme for the two-dimensional case. Section 5 provides some numerical results to clarify the schemes. The conclusion is given in the last section.

2. Legendre polynomials

The shifted Legendre polynomials $P^\sigma_i(x)$ , defined on the interval $\Omega = [0, \sigma]$ , are the eigenfunctions of the Sturm-Liouville problem:

$\frac{d}{dx} \left( x(\sigma - x) \frac{d}{dx} P^\sigma_i(x) \right) + i(i + 1) P^\sigma_i(x) = 0, \quad x \in \Omega,$

and can be expressed as:

$P^\sigma_i(x) = \sum\limits_{p = 0}^i C_{i, p} \frac{x^p}{\sigma^p} = \sum\limits_{p = 0}^i (-1)^{i-p} \frac{(i + p)!}{(i - p)! (p!)^2 \sigma^p} x^p,$

where $C_{i, p} = (-1)^{i-p} \frac{(i + p)!}{(i - p)! (p!)^2 \sigma^p}$ . These polynomials satisfy the derivative conditions:

$\frac{d^s}{dx^s} P^\sigma_i(0) = (-1)^{i-s} \frac{(i + s)!}{(i - s)! s! \sigma^s}, \quad \text{and} \quad \frac{d^s}{dx^s} P^\sigma_i(\sigma) = \frac{(i + s)!}{(i - s)! s! \sigma^s}.$

The shifted Legendre polynomials satisfy the orthogonality relation:

$\begin{equation} \int_0^\sigma P^\sigma_i(x) P^\sigma_p(x) \, dx = \frac{\sigma}{2i + 1} \delta_{ip}, \end{equation}$

(2.1)

where $\delta_{ip}$ is the Kronecker delta. Any function $y \in L^2(\Omega)$ can be expanded as:

$\begin{equation} y(x) = \sum\limits_{i = 0}^\infty a_i P^\sigma_i(x), \quad \text{where} \quad a_i = \frac{2i + 1}{\sigma} \int_0^\sigma y(x) P^\sigma_i(x) \, dx. \end{equation}$

(2.2)

Let $\text{Proj}^\sigma_{\mathcal{N}}$ denote the orthogonal projection operator:

$\text{Proj}^\sigma_{\mathcal{N}}: L^2(\Omega) \rightarrow S_{\mathcal{N}}, \quad S_{\mathcal{N}} = \text{Span}\{ P^\sigma_i(x) : 0 \leq i \leq \mathcal{N} \}.$

This projection is given by:

$\begin{equation} \text{Proj}^\sigma_{\mathcal{N}} y = y_{\mathcal{N}}(x) = \sum\limits_{i = 0}^{\mathcal{N}} a_i P^\sigma_i(x) = A_{\mathcal{N}}^T \Phi^\sigma_{\mathcal{N}}(x), \end{equation}$

(2.3)

where:

$\begin{equation} A_{\mathcal{N}} = \left[ a_i \quad \text{for} \quad 0 \leq i \leq \mathcal{N} \right]^T, \quad \Phi^\sigma_{\mathcal{N}}(x) = \left[ P^\sigma_i(x) \quad \text{for} \quad 0 \leq i \leq \mathcal{N} \right]^T. \end{equation}$

(2.4)

For functions $y \in L^2(\Omega_1 \times \Omega_2)$ with $\Omega_1 = [0, \sigma_1]$ and $\Omega_2 = [0, \sigma_2]$ , the approximation becomes:

$\begin{equation} y_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2) = \sum\limits_{i = 0}^{\mathcal{N}_1} \sum\limits_{j = 0}^{\mathcal{N}_2} a_{i, j} P^{\sigma_1}_i(x_1) P^{\sigma_2}_j(x_2) = A_{\mathcal{N}_1, \mathcal{N}_2}^T \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2), \end{equation}$

(2.5)

where:

$\begin{equation} \begin{split} A_{\mathcal{N}_1, \mathcal{N}_2} & = \left[ a_{i, j} \quad \text{for} \quad 0 \leq i \leq \mathcal{N}_1, \, 0 \leq j \leq \mathcal{N}_2 \right]^T, \\ \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2) & = \left[ P^{\sigma_1}_i(x_1) P^{\sigma_2}_j(x_2) \quad \text{for} \quad 0 \leq i \leq \mathcal{N}_1, \, 0 \leq j \leq \mathcal{N}_2 \right]^T. \end{split} \end{equation}$

(2.6)

The coefficients $a_{i, j}$ are computed via:

$\begin{equation} a_{i, j} = \frac{(2i + 1)(2j + 1)}{\sigma_1 \sigma_2} \int_0^{\sigma_1} \int_0^{\sigma_2} y(x_1, x_2) P^{\sigma_1}_i(x_1) P^{\sigma_2}_j(x_2) \, dx_1 dx_2. \end{equation}$

(2.7)

3. One-dimensional case

The spectral tau approach for (1.1) is to find $y_{\mathcal{N}}\in P_{\mathcal{N}}$ , such that

$\begin{equation} \begin{split} \sum\limits_{k = 0}^{n}b_k y_{\mathcal{N}}^{(k)}(x) = g_{\mathcal{N}}(x)+\sum\limits_{k = 0}^{n}c_k y_{\mathcal{N}}^{(k)}(p_k x)+\sum\limits_{k = 0}^{n}\int_{0}^{p_k x}d_k y_{\mathcal{N}}^{(k)}(s)ds+\int_{0}^{x}e y_{\mathcal{N}}(s)ds+\int_{0}^{\sigma}f y_{\mathcal{N}}(s)ds. \end{split} \end{equation}$

(3.1)

Now, we denote

$\begin{equation} \begin{split} y_{\mathcal{N}}(x)& = A^T_\mathcal{N} \Phi^{\sigma}_{\mathcal{N}}(x), \\ g_{\mathcal{N}}(x)& = \mathcal{G}^T_\mathcal{N} \Phi^{\sigma}_{\mathcal{N}}(x), \\ \end{split} \end{equation}$

(3.2)

where $\mathcal{G}_\mathcal{N}$ is a known vector given by

$\mathcal{G}_\mathcal{N} = \left[g_0, g_1, \cdots, g_\mathcal{N}\right]^T; \qquad g_j = \frac{2j+1}{\sigma} \int_{0}^{\sigma} g(x) P^{\sigma}_{j}(x) dx.$

If we denote $x_{\mathcal{N}, j, \sigma}$ by the Legendre Gauss quadrature nodes on the $(0, \sigma)$ and $\varpi_{\mathcal{N}, j, \sigma}$ by its corresponding Christoffel numbers, then we have

$g_j\simeq \frac{2j+1}{\sigma} \sum\limits_{k = 0}^{\mathcal{N}} \varpi_{\mathcal{N}, k, \sigma}g(x_{\mathcal{N}, k, \sigma}) P^{\sigma}_{j}(x_{\mathcal{N}, k, \sigma}).$

The following theorems are of great benefit later.

Theorem 1. ^[23] The derivative of order $k$ for the vector $\Phi^{\sigma}_{\mathcal{N}}(x)$ is given by:

$\begin{equation} \frac{d^k}{d x^k}\Phi^{\sigma}_{\mathcal{N}}(x) = \boldsymbol{D}_{(\mathcal{N}, k)}\Phi^{\sigma}_{\mathcal{N}}(x); \qquad \boldsymbol{D}_{(\mathcal{N}, k)} = \left(\boldsymbol{D}^{(1)}\right)^k, \end{equation}$

(3.3)

where $\boldsymbol{D}^{(1)} = (d_{r, j})_{0\leq r, j \leq \mathcal{N}},$ and

$d_{r, j} = \begin{cases}\dfrac{2}{\sigma}(2j+1), \ & r = j+s, \ \begin{cases}s = 1, 3, \cdots, \mathcal{N}-1, \ &\mathcal{N} \;\mathit{\text{is even}}, \\ s = 1, 3, \cdots, \mathcal{N}, \ &\mathcal{N} \;\mathit{\text{is odd}}, \end{cases}\\ 0, &otherwise. \end{cases}$

Theorem 2. ^[24] The integration of the vector $\Phi^{\sigma}_{\mathcal{N}}(x)$ is given by:

$\begin{equation} \int_0^{x} \Phi^{\sigma}_{\mathcal{N}}(x)dx = \boldsymbol{I}_{(\mathcal{N}, 1)} \Phi^{\sigma}_{\mathcal{N}}(x), \end{equation}$

(3.4)

where

$\boldsymbol{I}_{(\mathcal{N}, 1)} = (\beta_{i, j})_{0\leq i, j\leq \mathcal{N}}; \qquad \quad \beta_{i, j} = \sum\limits_{k = 0}^i\sum\limits_{l = 0}^{j}\frac{(-1)^{i+j+k+l} (i+k)!(j+l)!\sigma}{(i-k)!k!(k+1)!(j-l)!(l!)^2(k+l+2)}.$

Theorem 3. For $0 \leq p \leq 1,$ the pantograph operational matrix $\boldsymbol{Q}_{\mathcal{N}, p}$ can be defined by

$\begin{equation} \Phi^{\sigma}_{\mathcal{N}}(px) = \boldsymbol{Q}_{\mathcal{N}, p} \Phi^{\sigma}_{\mathcal{N}}(x), \end{equation}$

(3.5)

where $\boldsymbol{Q}_{\mathcal{N}, p} = (q^p_{k, j})_{0\leq k, j\leq \mathcal{N}},$ and $q^p_{k, j} = \sum\limits_{r = 0}^{k} \dfrac{C_{k, r}p^rf_{r, j}}{\sigma^r}$ .

Proof. We start by expressing $P^{\sigma}_{k}(px)$ by:

$\begin{equation} P^{\sigma}_{k}(px) = \sum\limits_{r = 0}^{k} C_{k, r} \frac{p^rx^r}{\sigma^r}. \end{equation}$

(3.6)

Expanding $x^r$ in terms of $P^{\sigma}_{j}(x), \ j = 0, 1, \cdots, \mathcal{N},$ by

$\begin{equation} x^r = \sum\limits_{j = 0}^{\mathcal{N}}f_{r, j}P^{\sigma}_{j}; \qquad \qquad f_{r, j} = \frac{2j+1}{\sigma} \int_{0}^{\sigma} x^r P^{\sigma}_{j}(x) dx, \qquad 0\leq j. \end{equation}$

(3.7)

A combination of (3.6) and (3.7) then yields

$\begin{equation} \begin{split} P^{\sigma}_{k}(px)& = \sum\limits_{r = 0}^{k} \frac{C_{k, r}p^r}{\sigma^r}\left(\sum\limits_{j = 0}^{\mathcal{N}}f_{r, j}P^{\sigma}_{j}\right) = \sum\limits_{j = 0}^{\mathcal{N}}P^{\sigma}_{j}\left(\sum\limits_{r = 0}^{k} \frac{C_{k, r}p^rf_{r, j}}{\sigma^r}\right)\\ & = \left[\sum\limits_{r = 0}^{k} \frac{C_{k, r}p^rf_{r, 0}}{\sigma^r}, \sum\limits_{r = 0}^{k} \frac{C_{k, r}p^rf_{r, 1}}{\sigma^r}, \cdots, \sum\limits_{r = 0}^{k} \frac{C_{k, r}p^rf_{r, \mathcal{N}}}{\sigma^r} \right]^T \Phi^{\sigma}_{\mathcal{N}}(x), \end{split} \end{equation}$

(3.8)

which completes the proof. □

Application of Theorems 1–3, we have

$\begin{equation} \begin{split} \frac{d^k y_{\mathcal{N}}(x)}{dx^k} & = A^T_\mathcal{N} \textbf{D}_{(\mathcal{N}, k)} \Phi^{\sigma}_{\mathcal{N}}(x), \\ \int_{0}^{x} y_{\mathcal{N}}(s)ds & = A^T_\mathcal{N} \textbf{I}_{(\mathcal{N}, 1)} \Phi^{\sigma}_{\mathcal{N}}(x), \\ y_{\mathcal{N}}(px)& = A^T_{\mathcal{N}}\textbf{Q}_{\mathcal{N}, p} \Phi^{\sigma}_{\mathcal{N}}(x), \\ \frac{d^k y_{\mathcal{N}}(p x)}{dx^k}& = A^T_{\mathcal{N}}\textbf{D}_{(\mathcal{N}, k)}\textbf{Q}_{\mathcal{N}, p} \Phi^{\sigma}_{\mathcal{N}}(x), \\ \int_{0}^{x} y_{\mathcal{N}}(p s)ds & = A^T_\mathcal{N} \textbf{I}_{(\mathcal{N}, 1)} \textbf{Q}_{\mathcal{N}, p} \Phi^{\sigma}_{\mathcal{N}}(x). \end{split} \end{equation}$

(3.9)

Using Theorems 1–3, and with the help of (3.9), the residual $\mathcal{R}_{\mathcal{N}}(x)$ for (3.1) is given by:

$\begin{equation} \begin{split} \mathcal{R}_{\mathcal{N}}(x)& = \sum\limits_{k = 0}^{n}b_k A^T_\mathcal{N} \textbf{D}_{(\mathcal{N}, k)} \Phi^{\sigma}_{\mathcal{N}}(x)-\sum\limits_{k = 0}^{n}c_k A^T_\mathcal{N} \textbf{D}_{(\mathcal{N}, k)} \textbf{Q}_{\mathcal{N}, p_k} \Phi^{\sigma}_{\mathcal{N}}(x)-e A^T_\mathcal{N} \textbf{I}_{(\mathcal{N}, 1)} \Phi^{\sigma}_{\mathcal{N}}(x)\\ &-\sum\limits_{k = 0}^{n}d_k A^T_\mathcal{N} \textbf{D}_{(\mathcal{N}, k)} \textbf{I}_{(\mathcal{N}, 1)} \textbf{Q}_{\mathcal{N}, p_k} \Phi^{\sigma}_{\mathcal{N}}(x)-f a_0\sigma-\mathcal{G}^T_\mathcal{N} \Phi^{\sigma}_{\mathcal{N}}(x), \end{split} \end{equation}$

(3.10)

According to the spectral tau method, the numerical solution for (1.1) is obtained by generating $\mathcal{N}+1$ linear algebraic equations as follows:

$\begin{equation} \begin{split} \int_{0}^{\sigma} \mathcal{R}_{\mathcal{N}}(x) P^{\sigma}_{k}(x) dx& = 0, \qquad \qquad k = 0, 1, \cdots, \mathcal{N}-n, \\ \end{split} \end{equation}$

(3.11)

$\begin{equation} \begin{split} A^T_\mathcal{N} \textbf{D}_{(\mathcal{N}, k)} \Phi^{\sigma}_{\mathcal{N}}(0)& = \varrho_k, \qquad \qquad k = 0, 1, \cdots, n-1. \end{split} \end{equation}$

(3.12)

Define the vector ${}_{\sigma}\textbf{M}_{\mathcal{N}}^i, \ i = 0, 1, \cdots, \mathcal{N}-n,$ as follows:

${}_{\sigma}\textbf{M}_{\mathcal{N}}^i = \frac{\sigma}{2i+1}\mathbf{e}_i,$

where $\mathbf{e}_i$ is the $\mathcal{N}+1$ standard basis vector. Then the system of equations (3.11) simplifies to:

$\begin{equation} \begin{split} \sum\limits_{k = 0}^{n}b_k A^T_\mathcal{N} \textbf{D}_{(\mathcal{N}, k)} & {}_{\sigma}\textbf{M}_{\mathcal{N}}^i-\sum\limits_{k = 0}^{n}c_k A^T_\mathcal{N} \textbf{D}_{(\mathcal{N}, k)} \textbf{Q}_{\mathcal{N}, p_k} {}_{\sigma}\textbf{M}_{\mathcal{N}}^i -\sum\limits_{k = 0}^{n}d_k A^T_\mathcal{N} \textbf{D}_{(\mathcal{N}, k)} \textbf{I}_{(\mathcal{N}, 1)} \textbf{Q}_{\mathcal{N}, p_k} {}_{\sigma}\textbf{M}_{\mathcal{N}}^i\\ & -e A^T_\mathcal{N} \textbf{I}_{(\mathcal{N}, 1)} {}_{\sigma}\textbf{M}_{\mathcal{N}}^i-f \sigma^2 \delta_{i0} A^T_\mathcal{N} \mathbf{e}_0 -\mathcal{G}^T_\mathcal{N} {}_{\sigma}\textbf{M}_{\mathcal{N}}^i = 0, \qquad i = 0, 1, \cdots, \mathcal{N}-n, \end{split} \end{equation}$

(3.13)

If we denote $E_i, \ i = 0, 1, \cdots, \mathcal{N}-n,$ by

$\begin{split} E_i = \sum\limits_{k = 0}^{n}b_k \left({}_{\sigma}\textbf{M}_{\mathcal{N}}^i\right)^T \textbf{D}_{(\mathcal{N}, k)}^T & -\sum\limits_{k = 0}^{n} c_k \left({}_{\sigma}\textbf{M}_{\mathcal{N}}^i\right)^T \textbf{Q}_{\mathcal{N}, p_k}^T \textbf{D}_{(\mathcal{N}, k)}^T -\sum\limits_{k = 0}^{n}d_k \left({}_{\sigma}\textbf{M}_{\mathcal{N}}^i\right)^T \textbf{Q}_{\mathcal{N}, p_k}^T \textbf{I}_{(\mathcal{N}, 1)}^T \textbf{D}_{(\mathcal{N}, k)}^T \\ & -e \left({}_{\sigma}\textbf{M}_{\mathcal{N}}^i\right)^T \textbf{I}_{(\mathcal{N}, 1)}^T -f \delta_{i0} \sigma^2 \mathbf{e}^T_0, \qquad i = 0, 1, \cdots, \mathcal{N}-n, \end{split}$

and

$C_i = \left(\Phi^{\sigma}_{\mathcal{N}}(0)\right)^T\textbf{D}_{(\mathcal{N}, i)}^T, \qquad i = 0, 1, \cdots, n-1,$

then the solution of the main problem is reduced to the system

$\textbf{S} A_{\mathcal{N}} = \textbf{B}.$

The coefficient matrix S is expressed as:

$\textbf{S} = \begin{bmatrix} \left({}_{\sigma}\textbf{M}_{\mathcal{N}}^0\right)^T \left( \sum\nolimits_{k = 0}^{n} b_k \textbf{D}_{(\mathcal{N}, k)}^T - c_k \textbf{Q}_{\mathcal{N}, p_k}^T \textbf{D}_{(\mathcal{N}, k)}^T - d_k \textbf{Q}_{\mathcal{N}, p_k}^T \textbf{I}_{(\mathcal{N}, 1)}^T \textbf{D}_{(\mathcal{N}, k)}^T - e \textbf{I}_{(\mathcal{N}, 1)}^T \right) -f \sigma^2 \mathbf{e}_0^T \\ \left({}_{\sigma}\textbf{M}_{\mathcal{N}}^1\right)^T \left( \sum\nolimits_{k = 0}^{n} b_k \textbf{D}_{(\mathcal{N}, k)}^T - c_k \textbf{Q}_{\mathcal{N}, p_k}^T \textbf{D}_{(\mathcal{N}, k)}^T - d_k \textbf{Q}_{\mathcal{N}, p_k}^T \textbf{I}_{(\mathcal{N}, 1)}^T \textbf{D}_{(\mathcal{N}, k)}^T - e \textbf{I}_{(\mathcal{N}, 1)}^T \right) \\ \vdots \\ \left({}_{\sigma}\textbf{M}_{\mathcal{N}}^{\mathcal{N}-n}\right)^T \left( \sum\nolimits_{k = 0}^{n} b_k \textbf{D}_{(\mathcal{N}, k)}^T - c_k \textbf{Q}_{\mathcal{N}, p_k}^T \textbf{D}_{(\mathcal{N}, k)}^T - d_k \textbf{Q}_{\mathcal{N}, p_k}^T \textbf{I}_{(\mathcal{N}, 1)}^T \textbf{D}_{(\mathcal{N}, k)}^T - e \textbf{I}_{(\mathcal{N}, 1)}^T \right) \\ \left( \Phi^{\sigma}_{\mathcal{N}}(0) \right)^T \textbf{D}_{(\mathcal{N}, 0)}^T \\ \vdots \\ \left( \Phi^{\sigma}_{\mathcal{N}}(0) \right)^T \textbf{D}_{(\mathcal{N}, n-1)}^T \end{bmatrix}.$

Compactly, we can write S as:

$\textbf{S} = \begin{bmatrix} E_0 \\ E_1 \\ \vdots \\ E_{\mathcal{N}-n} \\ C_0 \\ C_1 \\ \vdots \\ C_{n-1} \end{bmatrix},$

where each $E_i$ is a row vector of length $\mathcal{N}+1$ , and the total system is of size $(\mathcal{N}+1)\times (\mathcal{N}+1)$ . The right-hand side vector B is given by:

$\textbf{B} = \begin{bmatrix} \frac{g_0 \sigma}{1} \\ \frac{g_1 \sigma}{3} \\ \frac{g_2 \sigma}{5} \\ \vdots \\ \frac{g_{\mathcal{N}-n} \sigma}{2\mathcal{N}-2n+1} \\ \varrho_0 \\ \varrho_1 \\ \vdots \\ \varrho_{n-1} \end{bmatrix}.$

4. Two-dimensional case

In this section, we employ the Legendre spectral tau method to solve the two-dimensional high-order pantograph-type Volterra-Fredholm integro-differential equation:

$\begin{equation} \begin{split} &\sum\limits_{k_1 = 0}^{n_1}\sum\limits_{k_2 = 0}^{n_2} b_{k_1, k_2} \frac{\partial^{k_1+k_2} y(x_1, x_2)}{\partial x_1^{k_1} \partial x_2^{k_2}} = g(x_1, x_2) + \int_{0}^{\sigma_2} \int_{0}^{\sigma_1} c \, y(s_1, s_2) \, ds_1 ds_2 \\ &+ \int_{0}^{x_2} \int_{0}^{x_1} d \, y(s_1, s_2) \, ds_1 ds_2 + \sum\limits_{k_1 = 0}^{n_1} \sum\limits_{k_2 = 0}^{n_2} \int_{0}^{p_1 x_1} \int_{0}^{p_2 x_2} e_{k_1, k_2} \frac{\partial^{k_1+k_2} y(s_1, s_2)}{\partial x_1^{k_1} \partial x_2^{k_2}} \, ds_1 ds_2 \\ &+ \sum\limits_{k_1 = 0}^{n_1} \sum\limits_{k_2 = 0}^{n_2} f_{k_1, k_2} \frac{\partial^{k_1+k_2} y(p_1 x_1, p_2 x_2)}{\partial x_1^{k_1} \partial x_2^{k_2}}, \end{split} \end{equation}$

(4.1)

subject to the initial conditions:

$\begin{equation} \begin{split} \frac{\partial^{k_1} y(x_1, x_2)}{\partial x_1^{k_1}} & = \varrho_{1, k_1}(0, x_2), \quad k_1 = 0, 1, \dots, n_1-1, \quad 0 \leq x_2 \leq \sigma_2, \\ \frac{\partial^{k_2} y(x_1, x_2)}{\partial x_2^{k_2}} & = \varrho_{2, k_2}(x_1, 0), \quad k_2 = 0, 1, \dots, n_2-1, \quad 0 \leq x_1 \leq \sigma_1, \end{split} \end{equation}$

(4.2)

where $\sigma_1, \sigma_2, c, d, b_{k_1, k_2}, e_{k_1, k_2}$ , and $f_{k_1, k_2}$ (for $0 \leq k_1 \leq n_1$ , $0 \leq k_2 \leq n_2$ ) are known real constants.

As a spectral approach, we seek an approximate solution $y_{\mathcal{N}_1, \mathcal{N}_2} \in P_{\mathcal{N}_1} \times P_{\mathcal{N}_2}$ such that:

$\begin{equation} \begin{split} &\sum\limits_{k_1 = 0}^{n_1} \sum\limits_{k_2 = 0}^{n_2} b_{k_1, k_2} \frac{\partial^{k_1+k_2} y_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2)}{\partial x_1^{k_1} \partial x_2^{k_2}} = g_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2) \\ &+ \int_{0}^{x_2} \int_{0}^{x_1} d \, y_{\mathcal{N}_1, \mathcal{N}_2}(s_1, s_2) \, ds_1 ds_2 + \int_{0}^{\sigma_2} \int_{0}^{\sigma_1} c \, y_{\mathcal{N}_1, \mathcal{N}_2}(s_1, s_2) \, ds_1 ds_2 \\ &+ \sum\limits_{k_1 = 0}^{n_1} \sum\limits_{k_2 = 0}^{n_2} f_{k_1, k_2} \frac{\partial^{k_1+k_2} y_{\mathcal{N}_1, \mathcal{N}_2}(p_1 x_1, p_2 x_2)}{\partial x_1^{k_1} \partial x_2^{k_2}} \\ &+ \sum\limits_{k_1 = 0}^{n_1} \sum\limits_{k_2 = 0}^{n_2} \int_{0}^{p_1 x_1} \int_{0}^{p_2 x_2} e_{k_1, k_2} \frac{\partial^{k_1+k_2} y_{\mathcal{N}_1, \mathcal{N}_2}(s_1, s_2)}{\partial x_1^{k_1} \partial x_2^{k_2}} \, ds_1 ds_2. \end{split} \end{equation}$

(4.3)

We represent the approximate solution $y_{\mathcal{N}_1, \mathcal{N}_2}$ and the right-hand side function $g_{\mathcal{N}_1, \mathcal{N}_2}$ in terms of Legendre basis functions:

$\begin{equation} \begin{split} y_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2) & = A^T_{\mathcal{N}_1, \mathcal{N}_2} \, \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2), \\ g_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2) & = \mathcal{G}^T_{\mathcal{N}_1, \mathcal{N}_2} \, \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2), \end{split} \end{equation}$

(4.4)

where $\mathcal{G}_{\mathcal{N}_1, \mathcal{N}_2}$ is a known vector given by:

$\mathcal{G}_{\mathcal{N}_1, \mathcal{N}_2} = \left[ g_{0, 0}, g_{1, 0}, \dots, g_{\mathcal{N}_1, 0}, g_{0, 1}, \dots, g_{\mathcal{N}_1, \mathcal{N}_2} \right]^T,$

and the coefficients $g_{i, j}$ are computed as:

$g_{i, j} = \frac{(2i+1)(2j+1)}{\sigma_1 \sigma_2} \int_{0}^{\sigma_2} \int_{0}^{\sigma_1} g(x_1, x_2) \, P^{\sigma_1}_i(x_1) \, P^{\sigma_2}_j(x_2) \, dx_1 dx_2.$

The following matrix operators are defined for differentiation, integration, and scaling:

Theorem 4. Let $\mathbb{I}_{\mathcal{N}_1}$ and $\mathbb{I}_{\mathcal{N}_2}$ denote the identity matrices of orders ( $\mathcal{N}_1$ +1) and $(\mathcal{N}_2+1)$ , respectively. Then:

$\begin{equation} \frac{\partial^k}{\partial x_1^k} \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2) = \mathcal{D}_{(1, k)} \, \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2), \quad \mathcal{D}_{(1, k)} = \mathbf{D}_{(\mathcal{N}_1, 1)} \otimes \mathbb{I}_{\mathcal{N}_2}, \end{equation}$

(4.5)

$\begin{equation} \frac{\partial^k}{\partial x_2^k} \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2) = \mathcal{D}_{(2, k)} \, \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2), \quad \mathcal{D}_{(2, k)} = \mathbb{I}_{\mathcal{N}_1} \otimes \mathbf{D}_{(\mathcal{N}_2, 2)}, \end{equation}$

(4.6)

$\begin{equation} \int_{0}^{x_1} \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(s_1, x_2) \, ds_1 = \mathcal{I}_{(1, 1)} \, \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2), \quad \mathcal{I}_{(1, 1)} = \mathbf{I}_{(\mathcal{N}_1, 1)} \otimes \mathbb{I}_{\mathcal{N}_2}, \end{equation}$

(4.7)

$\begin{equation} \int_{0}^{x_2} \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, s_2) \, ds_2 = \mathcal{I}_{(2, 1)} \, \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2), \quad \mathcal{I}_{(2, 1)} = \mathbb{I}_{\mathcal{N}_1} \otimes \mathbf{I}_{(\mathcal{N}_2, 2)}, \end{equation}$

(4.8)

$\begin{equation} \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(p_1 x_1, x_2) = \mathcal{Q}_{1, p_1} \, \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2), \quad \mathcal{Q}_{1, p_1} = \mathbf{Q}_{\mathcal{N}_1, p_1} \otimes \mathbb{I}_{\mathcal{N}_2}, \end{equation}$

(4.9)

$\begin{equation} \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, p_2 x_2) = \mathcal{Q}_{2, p_2} \, \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2), \quad \mathcal{Q}_{2, p_2} = \mathbb{I}_{\mathcal{N}_1} \otimes \mathbf{Q}_{\mathcal{N}_2, p_2}. \end{equation}$

(4.10)

In virtue of (4.4)–(4.10), we have

$\begin{equation} \begin{split} \frac{\partial^k}{\partial x_1^k} y_{\mathcal{N}_1, \mathcal{N}_2} (x_1, x_2) & = A^T_{\mathcal{N}_1, \mathcal{N}_2} \mathcal{D}_{(1, k)}\Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2), \\ \frac{\partial^k}{\partial x_2^k} y_{\mathcal{N}_1, \mathcal{N}_2} (x_1, x_2) & = A^T_{\mathcal{N}_1, \mathcal{N}_2} \mathcal{D}_{(2, k)}\Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2), \\ \int_{0}^{x_1} y_{\mathcal{N}_1, \mathcal{N}_2} (s_1, x_2) ds_1 & = A^T_{\mathcal{N}_1, \mathcal{N}_2} \mathcal{I}_{(1, 1)}\Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2), \\ \int_{0}^{x_2} y_{\mathcal{N}_1, \mathcal{N}_2} (x_1, s_2) ds_2 & = A^T_{\mathcal{N}_1, \mathcal{N}_2} \mathcal{I}_{(2, 1)}\Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2), \\ y_{\mathcal{N}_1, \mathcal{N}_2} (p_1 x_1, x_2) & = A^T_{\mathcal{N}_1, \mathcal{N}_2} \mathcal{Q}_{1, p_1}\Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2), \\ y_{\mathcal{N}_1, \mathcal{N}_2} (x_1, p_2x_2) & = A^T_{\mathcal{N}_1, \mathcal{N}_2} \mathcal{Q}_{2, p_2}\Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2), \\ y_{\mathcal{N}_1, \mathcal{N}_2} (p_1x_1, p_2x_2) & = A^T_{\mathcal{N}_1, \mathcal{N}_2} \mathcal{Q}_{1, p_2}\mathcal{Q}_{2, p_2}\Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2). \\ \end{split} \end{equation}$

(4.11)

Using the basis representation and matrix operators, the residual $\mathcal{R}_{\mathcal{N}_1, \mathcal{N}_2}$ of (4.3) is given by:

$\begin{split} \mathcal{R}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2) & = \sum\limits_{k_1 = 0}^{n_1} \sum\limits_{k_2 = 0}^{n_2} b_{k_1, k_2} \, A^T_{\mathcal{N}_1, \mathcal{N}_2} \, \mathcal{D}_{(1, k_1)} \, \mathcal{D}_{(2, k_2)} \, \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2) - \mathcal{G}^T_{\mathcal{N}_1, \mathcal{N}_2} \, \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2)\\ &- d \, A^T_{\mathcal{N}_1, \mathcal{N}_2} \, \mathcal{I}_{(1, 1)} \, \mathcal{I}_{(2, 1)} \, \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2) - c \sigma_1\sigma_2a_{0, 0} \\ &- \sum\limits_{k_1 = 0}^{n_1} \sum\limits_{k_2 = 0}^{n_2} e_{k_1, k_2} \, A^T_{\mathcal{N}_1, \mathcal{N}_2} \, \mathcal{D}_{(1, k_1)} \, \mathcal{D}_{(2, k_2)} \, \mathcal{I}_{(1, 1)} \, \mathcal{I}_{(2, 1)} \, \mathcal{Q}_{1, p_1} \, \mathcal{Q}_{2, p_2} \, \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2) \\ &- \sum\limits_{k_1 = 0}^{n_1} \sum\limits_{k_2 = 0}^{n_2} f_{k_1, k_2} \, A^T_{\mathcal{N}_1, \mathcal{N}_2} \, \mathcal{D}_{(1, k_1)} \, \mathcal{D}_{(2, k_2)} \, \mathcal{Q}_{1, p_1} \, \mathcal{Q}_{2, p_2} \, \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2). \end{split}$

The system of $(\mathcal{N}_1+1)(\mathcal{N}_2+1)$ algebraic equations is obtained by enforcing the residual to be orthogonal to the basis functions:

$\begin{equation} \begin{split} &\int_{0}^{\sigma_2} \int_{0}^{\sigma_1} \mathcal{R}_{\mathcal{N}_1, \mathcal{N}_2}(x_{1, i}, x_{2, j}) \, P^{\sigma_1}_i(x_1) \, P^{\sigma_2}_j(x_2) \, dx_1 dx_2 = 0, \\ &\quad 0 \leq i \leq \mathcal{N}_1 - n_1, \quad 0 \leq j \leq \mathcal{N}_2 - n_2, \\ \end{split} \end{equation}$

(4.12)

$\begin{equation} \begin{split} &A^T_{\mathcal{N}_1, \mathcal{N}_2} \, \mathcal{D}_{(1, k_1)} \, \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(0, x_{2, j}) = \varrho_{1, k_1}(x_{2, j}), \quad 0 \leq k_1 \leq n_1-1, \quad 0 \leq j \leq \mathcal{N}_2, \\ \end{split} \end{equation}$

(4.13)

$\begin{equation} \begin{split} &A^T_{\mathcal{N}_1, \mathcal{N}_2} \, \mathcal{D}_{(2, k_2)} \, \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_{1, j}, 0) = \varrho_{2, k_2}(x_{1, j}), \quad 0 \leq k_2 \leq n_2-1, \quad n_1 \leq j \leq \mathcal{N}_1, \end{split} \end{equation}$

(4.14)

where $x_{1, i}$ and $x_{2, j}$ are the roots of $P^{\sigma_1}_{\mathcal{N}_1+1}(x_1)$ and $P^{\sigma_2}_{\mathcal{N}_2+1}(x_2)$ , respectively.

Define the vector $\mathbf{e}_{00} = \mathbf{e}_0\otimes\mathbf{e}_0$ and ${}_{\sigma_1, \sigma_2}\textbf{M}_{\mathcal{N}_1, \mathcal{N}_2}^{i, j}, \ i = 0, 1, \cdots, \mathcal{N}_1-n_1, \ j = 0, 1, \cdots, \mathcal{N}_2-n_2,$ as follows:

${}_{\sigma_1, \sigma_2} \textbf{M}_{\mathcal{N}_1, \mathcal{N}_2}^{i, j} = {}_{\sigma_1} \textbf{M}^i_{\mathcal{N}_1}\otimes {}_{\sigma_2} \textbf{M}^j_{\mathcal{N}_2}, \quad i = 0, 1, \cdots, \mathcal{N}_1-n_1, \ j = 0, 1, \cdots, \mathcal{N}_2-n_2,$

then the system (4.12) is simplified to:

$\begin{split} \sum\limits_{k_1 = 0}^{n_1} \sum\limits_{k_2 = 0}^{n_2} & b_{k_1, k_2} \, A^T_{\mathcal{N}_1, \mathcal{N}_2} \, \mathcal{D}_{(1, k_1)} \, \mathcal{D}_{(2, k_2)} \, {}_{\sigma_1, \sigma_2}\textbf{M}_{\mathcal{N}_1, \mathcal{N}_2}^{i, j} - \mathcal{G}^T_{\mathcal{N}_1, \mathcal{N}_2} \, {}_{\sigma_1, \sigma_2}\textbf{M}_{\mathcal{N}_1, \mathcal{N}_2}^{i, j}- d \, A^T_{\mathcal{N}_1, \mathcal{N}_2} \, \mathcal{I}_{(1, 1)} \, \mathcal{I}_{(2, 1)} \, {}_{\sigma_1, \sigma_2}\textbf{M}_{\mathcal{N}_1, \mathcal{N}_2}^{i, j}\\ &- c \sigma_1^2 \sigma_2^2 \delta_{i0} \delta_{j0} A^T_{\mathcal{N}_1, \mathcal{N}_2} \, \mathbf{e}_{00} - \sum\limits_{k_1 = 0}^{n_1} \sum\limits_{k_2 = 0}^{n_2} f_{k_1, k_2} \, A^T_{\mathcal{N}_1, \mathcal{N}_2} \, \mathcal{D}_{(1, k_1)} \, \mathcal{D}_{(2, k_2)} \, \mathcal{Q}_{1, p_1} \, \mathcal{Q}_{2, p_2} \, {}_{\sigma_1, \sigma_2}\textbf{M}_{\mathcal{N}_1, \mathcal{N}_2}^{i, j} \\ &- \sum\limits_{k_1 = 0}^{n_1} \sum\limits_{k_2 = 0}^{n_2} e_{k_1, k_2} \, A^T_{\mathcal{N}_1, \mathcal{N}_2} \, \mathcal{D}_{(1, k_1)} \, \mathcal{D}_{(2, k_2)} \, \mathcal{I}_{(1, 1)} \, \mathcal{I}_{(2, 1)} \, \mathcal{Q}_{1, p_1} \, \mathcal{Q}_{2, p_2} \, {}_{\sigma_1, \sigma_2}\textbf{M}_{\mathcal{N}_1, \mathcal{N}_2}^{i, j}. \end{split}$

Denoting $E_{i, j}, \ i = 0, 1, \cdots, \mathcal{N}_1-n_1, \ j = 0, 1, \cdots, \mathcal{N}_2-n_2,$ by

$\begin{split} E_{i, j}& = \sum\limits_{k_1 = 0}^{n_1} \sum\limits_{k_2 = 0}^{n_2} b_{k_1, k_2} \, \left({}_{\sigma_1, \sigma_2}\textbf{M}_{\mathcal{N}_1, \mathcal{N}_2}^{i, j}\right)^T \, \mathcal{D}^T_{(2, k_2)} \, \mathcal{D}^T_{(1, k_1)} - d \, \left({}_{\sigma_1, \sigma_2}\textbf{M}_{\mathcal{N}_1, \mathcal{N}_2}^{i, j}\right)^T \mathcal{I}^T_{(2, 1)} \mathcal{I}^T_{(1, 1)}\\ &- c \, \sigma_1^2 \sigma_2^2 \delta_{i0} \delta_{j0} \mathbf{e}_{00}^T - \sum\limits_{k_1 = 0}^{n_1} \sum\limits_{k_2 = 0}^{n_2} f_{k_1, k_2} \, \left({}_{\sigma_1, \sigma_2}\textbf{M}_{\mathcal{N}_1, \mathcal{N}_2}^{i, j}\right)^T \mathcal{Q}^T_{2, p_2} \mathcal{Q}^T_{1, p_1} \mathcal{D}^T_{(2, k_2)}\mathcal{D}^T_{(1, k_1)} \\ &- \sum\limits_{k_1 = 0}^{n_1} \sum\limits_{k_2 = 0}^{n_2} e_{k_1, k_2} \, \left({}_{\sigma_1, \sigma_2}\textbf{M}_{\mathcal{N}_1, \mathcal{N}_2}^{i, j}\right)^T \mathcal{Q}^T_{2, p_2} \mathcal{Q}^T_{1, p_1} \mathcal{I}^T_{(2, 1)} \mathcal{I}^T_{(1, 1)} \mathcal{D}^T_{(2, k_2)} \mathcal{D}^T_{(1, k_1)} , \end{split}$

and

$\begin{split} &C_{1, i, j} = \left(\Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(0, x_{2, j}) \right)^T\mathcal{D}^T_{(1, i)}, \qquad i = 0, 1, \cdots, n_1-1, \quad j = 0, 1, \cdots, \mathcal{N}_2, \\ &C_{2, i, j} = \left(\Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_{1, j}, 0)\right)^T\mathcal{D}^T_{(2, i)}, \qquad i = 0, 1, \cdots, n_2-1, \quad j = n_1, n_1+1, \cdots, \mathcal{N}_1. \end{split}$

The solution of the main problem is reduced to the system

$\textbf{S} A_{\mathcal{N}_1, \mathcal{N}_2} = \textbf{B},$

where

$\begin{split} \textbf{S} = \Big[E_{0, 0};\ E_{0, 1}&;\ \cdots;\ E_{\mathcal{N}_1-n_1;\mathcal{N}_2-n_2};\ C_{1, 0, 0};\ C_{1, 0, \mathcal{N}_2};\ \cdots;\ C_{1, n_1-1, \mathcal{N}_2};\ C_{2, 0, n_1};\ C_{2, 0, n_1+1};\ \cdots;\ C_{2, n_2-1, \mathcal{N}_1}\Big]^T, \\ \textbf{B} = \Big[\frac{g_{0, 0}\sigma_1\sigma_2}{(1)(1)}&, \ \frac{g_{0, 1}\sigma_1\sigma_2}{(1)(3)}, \frac{g_{0, 2}\sigma_1\sigma_2}{(1)(5)}, \ \cdots, \frac{g_{\mathcal{N}_1-n_1, \mathcal{N}_2-n_2}\sigma_1\gamma_2}{(2\mathcal{N}_1-2n_1+1)(2\mathcal{N}_2-2n_2+1)} , \varrho_{1, 0}(x_{2, 0}), \\ & \varrho_{1, 0}(x_{2, 1}), \ \cdots, \ \varrho_{1, n_1-1}(x_{2, \mathcal{N}_2}), \ \varrho_{2, 0}(x_{1, n_1}), \varrho_{2, 0}(x_{1, n_1+1}), \cdots, \ \varrho_{2, n_2-1}(x_{1, \mathcal{N}_1]}) \Big]^T. \end{split}$

5. Numerical results

5.1. Convergence test

In this subsection, we analyze the convergence behavior of the proposed numerical method. We consider a test problem with an irregular solution and compute the maximum absolute errors (MAEs) for different values of $\mathcal{N}$ .

We consider the following P-VF-IDE:

$\begin{equation} \frac{d^2 y(x)}{dx^2} + \frac{dy(x)}{dx} = \frac{1}{2} y\left(\frac{1}{3}x\right) + \int_{0}^{x} y(t) \, dt + \int_{0}^{\frac{1}{3}x} y(t) \, dt + \int_{0}^{1} y(t) \, dt + g(x), \quad 0 \leq x \leq 1, \end{equation}$

(5.1)

with the initial conditions:

$y(0) = 0, \quad \frac{dy(0)}{dx} = 0.$

The function $g(x)$ is chosen such that the exact solution is:

$y(x) = x^{\frac{11}{3}}.$

We apply the proposed numerical scheme to this problem for various values of $\mathcal{N}$ . The MAEs of the numerical solution $y_{\mathcal{N}}(x)$ are computed and summarized in . The MAEs decrease as $\mathcal{N}$ increases, demonstrating the method's high accuracy. The convergence behavior is further visualized in . The results demonstrate that the proposed numerical method achieves excellent convergence for problems with irregular solutions. The decay of errors with increasing $\mathcal{N}$ highlights the efficiency and robustness of the method.

Table 1. The MAEs of

$y_{\mathcal{N}}(x)$ for problem 5.1.

$\mathcal{N}$	4	8	12	16	20
MAE	$8.2705\times10^{-4}$	$4.9139\times10^{-6}$	$6.0563\times10^{-7}$	$1.5001\times10^{-7}$	$4.9759\times10^{-8}$

| Show Table

DownLoad: CSV

Figure 1. Convergence of the approximate solution for problem 5.1.

DownLoad: Full-Size Img PowerPoint

5.2. The numerical stability test

To numerically investigate the stability of the spectral tau method, we consider the following Volterra integro-differential equation ^[25]:

$\begin{equation} \begin{split} \frac{d y(x)}{dx} = y(0.5x) + \int_{0}^{x} y(t) \, dt + \int_{0}^{0.5x} y(t) \, dt + 1 - 1.5x, \quad 0 \leq x \leq 1, \end{split} \end{equation}$

(5.2)

with the initial condition $y(0) = 0$ and the exact solution $y(x) = 1 - e^{-x}$ .

To assess the stability of the method, we introduce perturbations to the right-hand side and the initial condition. Specifically, we consider the following perturbed problems:

$(i)$ Perturbed right-hand side: The perturbed right-hand side problem is given by:

$\begin{equation} \frac{d z(x)}{dx} = z(0.5x) + \int_{0}^{x} z(t) \, dt + \int_{0}^{0.5x} z(t) \, dt + 1 - 1.5x + \epsilon_r, \end{equation}$

(5.3)

with the initial condition $z(0) = 0$ and the exact solution $z(x) = 1 - e^{-x}$ , where $\epsilon_r$ is a small perturbation parameter. This perturbation tests the sensitivity of the method to changes in the forcing term of the equation.

The maximum absolute errors, $\lvert y_{\mathcal{N}}(x) - z_{\mathcal{N}}(x) \rvert$ , for the perturbed right-hand side problem (5.3) are computed for several values of $\epsilon_r$ . The results are summarized in Table 2.

Table 2. Maximum absolute errors for the perturbed right-hand side problem (5.3).

$\mathcal{N}$	$\epsilon_r=0.1$	$\epsilon_r=0.01$	$\epsilon_r=0.001$
6	0.15237708	0.01523770	0.00152376
8	0.15237709	0.01523770	0.00152377
10	0.15237709	0.01523770	0.00152377

| Show Table

DownLoad: CSV

$(ii)$ Perturbed initial condition: The perturbed initial condition problem is given by:

$\begin{equation} \frac{d z(x)}{dx} = z(0.5x) + \int_{0}^{x} z(t) \, dt + \int_{0}^{0.5x} z(t) \, dt + 1 - 1.5x, \end{equation}$

(5.4)

with the perturbed initial condition:

$z(0) = \epsilon_i,$

where $\epsilon_i$ is a small perturbation parameter. This perturbation tests the sensitivity of the method to changes in the initial condition.

The maximum absolute errors, $\lvert y_{\mathcal{N}}(x) - z_{\mathcal{N}}(x) \rvert$ , for the perturbed initial condition problem (5.4) are computed for several values of $\epsilon_i$ . The results are summarized in Table 3.

Table 3. Maximum absolute errors for the perturbed initial condition problem (5.4).

$\mathcal{N}$	$\epsilon_i=0.1$	$\epsilon_i=0.01$	$\epsilon_i =0.001$
6	0.34154212	0.03415420	0.00341541
8	0.34154212	0.03415421	0.00341542
10	0.34154212	0.03415421	0.00341542

| Show Table

DownLoad: CSV

The numerical results demonstrate that the spectral tau method remains stable under small perturbations to both the right-hand side and the initial condition. The errors introduced by the perturbations are proportional to the perturbation parameters $\epsilon_r$ and $\epsilon_i$ , indicating that the method is robust and well-conditioned. Specifically:

● For the perturbed right-hand side, the errors decrease linearly with $\epsilon_r$ .

● For the perturbed initial condition, the errors decrease linearly with $\epsilon_i$ .

These findings confirm the stability of the spectral tau method and its ability to handle small perturbations without significant loss of accuracy.

5.3. Comparison test

In this subsection, we evaluate the performance of the proposed numerical method by comparing it with existing methods. Specifically, we consider the following Volterra integro-differential equation ^[25]:

$\begin{equation} \begin{split} \frac{d y(x)}{dx} = y(0.5x) + \int_{0}^{x} y(t) \, dt + \int_{0}^{0.5x} y(t) \, dt + 1 - 1.5x, \quad 0 \leq x \leq 1, \end{split} \end{equation}$

(5.5)

with the initial condition $y(0) = 0$ and the exact solution $y(x) = 1 - e^{-x}$ .

For the solution of this problem, Zhao et al. ^[25] used the spectral collocation approach based on the Sinc function (SCSFM). This method reduces the problem to solving a system of algebraic equations. In , we compare the MAEs of $y_{\mathcal{N}}(x)$ obtained using our method with those obtained using the SCSFM ^[25].

Table 4. Comparing MAEs of

$y_{\mathcal{N}}(x)$ given using our method against those given using the SCSFM ^[25] for 5.5.

SCSFM ^[25]		Our scheme
$N$	MAE	$M$	MAE
10	$1.2608\times10^{-4}$	2	$7.1205\times10^{-3}$
20	$2.2598\times10^{-6}$	4	$1.6006\times10^{-5}$
30	$9.5984\times10^{-8}$	6	$2.2379\times10^{-8}$
40	$7.1028\times10^{-9}$	8	$1.9964\times10^{-11}$
50	$7.4603\times10^{-10}$	10	$1.0957\times10^{-14}$
60	$9.7609\times10^{-11}$	12	$4.0571\times10^{-16}$

| Show Table

DownLoad: CSV

The results in demonstrate that our method achieves significantly higher accuracy compared to the SCSFM, especially for larger values of $\mathcal{N}$ . For instance, with $\mathcal{N} = 12$ , our method achieves an MAE of $4.0571 \times 10^{-16}$ , whereas the SCSFM achieves $9.7609 \times 10^{-11}$ with $\mathcal{N} = 60$ . This highlights the superior convergence properties of our method.

To further illustrate the accuracy of our method, we plot the absolute error function $|y(x) - y_{\mathcal{N}}(x)|$ for $\mathcal{N} = 6, 10$ and $14$ in Figure 2.

Figure 2. Absolute errors with

$\mathcal{N} = 6, 10,$ and

$14$ for problem 5.5.

DownLoad: Full-Size Img PowerPoint

5.4. Two-dimensional case

In this subsection, we extend the application of the proposed numerical method to a two-dimensional problem. The goal is to demonstrate the effectiveness of the method in solving higher-dimensional integro-differential equations. We consider the following problem:

$\begin{equation} \begin{split} \frac{\partial^{5} y(x_1, x_2)}{\partial x_1^{3} \partial x_2^{2}} & = -y(x_1, x_2) + y(0.5x_1, 0.25x_2) - \int_{0}^{0.25x_2} \int_{0}^{0.5x_1} \frac{\partial y(s_1, s_2)}{\partial x_1} \, ds_1 ds_2 \\ &+ \int_{0}^{x_2} \int_{0}^{x_1} y(s_1, s_2) \, ds_1 ds_2 + \int_{0}^{1} \int_{0}^{1} y(s_1, s_2) \, ds_1 ds_2 + g(x_1, x_2), \end{split} \end{equation}$

(5.6)

where $0 \leq x_1, x_2 \leq 1$ , and the initial conditions are given by:

$y(0, x_2) = \log(x_2 + 1), \quad \frac{\partial y(0, x_2)}{\partial x_1} = \frac{\partial^2 y(0, x_2)}{\partial x_1^2} = y(x_1, 0) = 0, \quad \frac{\partial y(x_1, 0)}{\partial x_2} = x_1^4 + 1.$

The function $g(x_1, x_2)$ is chosen such that the exact solution is:

$y(x_1, x_2) = (x_1^4 + 1) \log(1 + x_2).$

We apply the numerical scheme presented in Section 4 to solve this problem. The discretization is performed with $\mathcal{N}_1 = 4$ and varying values of $\mathcal{N}_2$ . The absolute errors of the numerical solution $y_{4, \mathcal{N}_2}(x_1, x_2)$ are computed for $\mathcal{N}_2 = 4, 8, 12, 16,$ and $20$ . The results are summarized in Table 5.

Table 5. Absolute errors of

$y_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2)$ at

$\mathcal{N}_1 = 4$ and different choices of

$\mathcal{N}_2$ for problem 5.6.

$(x_1, x_2)$	$\mathcal{N}_2=4$	$\mathcal{N}_2=8$	$\mathcal{N}_2=12$	$\mathcal{N}_2=16$	$\mathcal{N}_2=20$
(0.1, 0.1)	$5.695\times10^{-5}$	$2.231\times10^{-8}$	$7.399\times10^{-12}$	$3.955\times10^{-15}$	$4.163\times10^{-17}$
(0.2, 0.2)	$1.763\times10^{-5}$	$7.956\times10^{-8}$	$1.358\times10^{-10}$	$1.518\times10^{-13}$	$3.330\times10^{-16}$
(0.3, 0.3)	$2.306\times10^{-4}$	$3.891\times10^{-7}$	$5.919\times10^{-10}$	$7.792\times10^{-13}$	$1.111\times10^{-15}$
(0.4, 0.4)	$6.320\times10^{-4}$	$1.224\times10^{-6}$	$1.783\times10^{-9}$	$2.277\times10^{-12}$	$2.997\times10^{-15}$
(0.5, 0.5)	$1.367\times10^{-3}$	$2.739\times10^{-6}$	$4.033\times10^{-9}$	$5.164\times10^{-12}$	$6.772\times10^{-15}$
(0.6, 0.6)	$2.603\times10^{-3}$	$5.225\times10^{-6}$	$7.698\times10^{-9}$	$9.859\times10^{-12}$	$1.276\times10^{-14}$
(0.7, 0.7)	$4.467\times10^{-3}$	$8.894\times10^{-6}$	$1.303\times10^{-8}$	$1.664\times10^{-11}$	$2.164\times10^{-14}$
(0.8, 0.8)	$7.000\times10^{-3}$	$1.369\times10^{-5}$	$2.004\times10^{-8}$	$2.563\times10^{-11}$	$3.330\times10^{-14}$
(0.9, 0.9)	$1.007\times10^{-2}$	$1.953\times10^{-5}$	$2.860\times10^{-9}$	$3.654\times10^{-11}$	$4.707\times10^{-14}$

| Show Table

DownLoad: CSV

The results in demonstrate the high accuracy of the proposed method for two-dimensional problems. As $\mathcal{N}_2$ increases, the absolute errors decrease significantly, indicating the convergence of the numerical solution to the exact solution. To further illustrate the accuracy of the method, we plot the absolute error function $|y(x_1, x_2) - y_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2)|$ for $(\mathcal{N}_1, \mathcal{N}_2) = \{(4, 10), (4, 14), (4, 18)\}$ in . To analyze the convergence behavior, we plot the logarithmic function of the MAEs for $\mathcal{N}_1 = 4$ and varying $\mathcal{N}_2$ in Figure 4.

Figure 3. Absolute error function of

$y(x_1, x_2)$ with

$(\mathcal{N}_1, \mathcal{N}_2) = \{(4, 10), (4, 14), (4, 18)\}$ for problem 5.6.

DownLoad: Full-Size Img PowerPoint

Figure 4. The logarithmic function of MAEs of

$y(x_1, x_2)$ with

$\mathcal{N}_1 = 4$ for problem 5.6.

DownLoad: Full-Size Img PowerPoint

The numerical results confirm the effectiveness of the proposed method for solving two-dimensional integro-differential equations. The errors decrease exponentially as $\mathcal{N}_2$ increases, demonstrating the spectral accuracy of the method. The logarithmic error plot further highlights the rapid convergence of the numerical solution.

6. Conclusions

This paper investigated a class of high-order P-VF-IDEs, which incorporate both Volterra and Fredholm integral components along with pantograph delay elements. We introduced a spectral tau approach for approximating solutions to P-VF-IDEs in one and two dimensions, utilizing operational differentiation and integration matrices to transform the continuous problem into a manageable system of algebraic equations. This method demonstrated high accuracy with a few numbers of computational modes. Through comprehensive numerical experiments, we highlighted the accuracy and convergence properties of the spectral Legendre tau method, affirming its effectiveness in solving high-order P-VF-IDEs in comparison to other spectral techniques. The results indicate that the proposed approach is a powerful tool for addressing intricate integro-differential equations with integral and proportional delay features.

Author contributions

M. A. Zaky: Writing–review & editing, Writing–original draft, Validation, Supervision, Software, Investigation; W. G. Alharbi: Validation, Formal analysis; M. M. Alzubaidi: Validation, Methodology, Writing the original draft. R. T. Matoog: Validation, Methodology, Writing the original draft. All authors have read and agreed to the published version of the manuscript.

Use of Generative-AI tools declaration

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

Acknowledgments

This work was supported and funded by the Deanship of Scientific Research at Imam Mohammad Ibn Saud Islamic University (IMSIU) (grant number IMSIU-DDRSP2503).

Funding

This work was supported and funded by the Deanship of Scientific Research at Imam Mohammad Ibn Saud Islamic University (IMSIU) (grant number IMSIU-DDRSP2503).

Conflict of interest

The authors declare that they do not have a conflict of interest.

References

[1]	R. A. Meyers, Encyclopedia of physical science and technology, San Diego: Academic Press, 1992.
[2]	C. Li, W. Hu, J. C. Sprott, X. Wang, Multistability in symmetric chaotic systems, Eur. Phys. J. Spec. Top., 224 (2015), 1493–1506. https://doi.org/10.1140/epjst/e2015-02475-x doi: 10.1140/epjst/e2015-02475-x
[3]	T. A. Alexeeva, N. V. Kuznetsov, T. N. Mokaev, Study of irregular dynamics in an economic model: attractor localization and Lyapunov exponents, Chaos Soliton. Fract., 152 (2021), 111365. https://doi.org/10.1016/j.chaos.2021.111365 doi: 10.1016/j.chaos.2021.111365
[4]	S. Vaidyanathan, A. S. T. Kammogne, E. Tlelo-Cuautle, C. N. Talonang, B. Abd-El-Atty, A. A. Abd El-Latif, et al., A novel 3-D jerk system, its bifurcation analysis, electronic circuit design and a cryptographic application, Electronics, 12 (2023), 2818. https://doi.org/10.3390/electronics12132818 doi: 10.3390/electronics12132818
[5]	C. Nwachioma, J. H. P'erez-Cruz, Analysis of a new chaotic system, electronic realization and use in navigation of differential drive mobile robot, Chaos Soliton. Fract., 144 (2021), 110684. https://doi.org/10.1016/j.chaos.2021.110684 doi: 10.1016/j.chaos.2021.110684
[6]	N. V. Kuznetsov, G. A. Leonov, V. I. Vagaitsev, Analytical-numerical method for attractor localization of generalized Chua's system, IFAC Proc. Vol., 43 (2010), 29–33. https://doi.org/10.3182/20100826-3-TR-4016.00009 doi: 10.3182/20100826-3-TR-4016.00009
[7]	G. A. Leonov, N. V. Kuznetsov, V. I. Vagaitsev, Localization of hidden Chua's attractors, Phys. Lett. A, 375 (2011), 2230–2233. https://doi.org/10.1016/j.physleta.2011.04.037 doi: 10.1016/j.physleta.2011.04.037
[8]	G. A. Leonov, N. V. Kuznetsov, Hidden attractors in dynamical systems. From hidden oscillations in Hilbert-Kolmogorov, Aizerman, and Kalman problems to hidden chaotic attractor in Chua circuits, Int. J. Bifurcat. Chaos, 23 (2013), 1330002. https://dx.doi.org/10.1142/S0218127413300024 doi: 10.1142/S0218127413300024
[9]	N. V. Stankevich, N. V. Kuznetsov, G. A. Leonov, L. O. Chua, Scenario of the birth of hidden attractors in the Chua circuit, Int. J. Bifurcat. Chaos, 27 (2017), 1730038. https://doi.org/10.1142/S0218127417300385 doi: 10.1142/S0218127417300385
[10]	N. Kuznetsov, T. Mokaev, V. Ponomarenko, E. Seleznev, N. Stankevich, L. Chua, Hidden attractors in Chua circuit: mathematical theory meets physical experiments, Nonlinear Dyn., 111 (2023), 5859–5887. https://doi.org/10.1007/s11071-022-08078-y doi: 10.1007/s11071-022-08078-y
[11]	Q. Wu, Q. Hong, X. Liu, X. Wang, Z. Zeng, A novel amplitude control method for constructing nested hidden multi-butterfly and multiscroll chaotic attractors, Chaos Soliton. Fract., 134 (2020), 109727. https://doi.org/10.1016/j.chaos.2020.109727 doi: 10.1016/j.chaos.2020.109727
[12]	H. Tian, Z. Wang, P. Zhang, M. Chen, Y. Wang, Dynamic analysis and robust control of a chaotic system with hidden attractor, Complexity, 2021 (2021), 8865522. https://doi.org/10.1155/2021/8865522 doi: 10.1155/2021/8865522
[13]	F. Bao, S. Yu, How to generate chaos from switching system: a saddle focus of index 1 and heteroclinic loop-based approach, Math. Probl. Eng., 2011 (2011), 756462. https://doi.org/10.1155/2011/756462 doi: 10.1155/2011/756462
[14]	T. Zhou, G. Chen, Classification of chaos in 3-D autonomous quadratic systems-Ⅰ: basic framework and methods, Int. J. Bifurcat. Chaos, 16 (2006), 2459–2479. https://doi.org/10.1142/S0218127406016203 doi: 10.1142/S0218127406016203
[15]	Z. Wei, J. C. Sprott, H. Chen, Elementary quadratic chaotic flows with a single non-hyperbolic equilibrium, Phys. Lett. A, 379 (2015), 2184–2187. https://doi.org/10.1016/j.physleta.2015.06.040 doi: 10.1016/j.physleta.2015.06.040
[16]	K. Rajagopal, V. T. Pham, F. R. Tahir, A. Akgul, H. R. Abdolmohammadi, S. Jafari, A chaotic jerk system with non-hyperbolic equilibrium: dynamics, effect of time delay and circuit realisation, Pramana, 90 (2018), 52. https://doi.org/10.1007/s12043-018-1545-x doi: 10.1007/s12043-018-1545-x
[17]	X. Cai, L. Liu, Y. Wang, C. Liu, A 3D chaotic system with piece-wise lines shape non-hyperbolic equilibria and its predefined-time control, Chaos Soliton. Fract., 146 (2021), 110904. https://doi.org/10.1016/j.chaos.2021.110904 doi: 10.1016/j.chaos.2021.110904
[18]	C. Li, W. Hai, Constructing multiwing attractors from a robust chaotic system with non-hyperbolic equilibrium points, Automatika, 59 (2018), 184–193. https://doi.org/10.1080/00051144.2018.1516273 doi: 10.1080/00051144.2018.1516273
[19]	C. Li, J. C. Sprott, Chaotic flows with a single nonquadratic term, Phys. Lett. A, 378 (2014), 178–183. https://doi.org/10.1016/j.physleta.2013.11.004 doi: 10.1016/j.physleta.2013.11.004
[20]	S. Jafari, J. C. Sprott, F. Nazarimehr, Recent new examples of hidden attractors, Eur. Phys. J. Spec. Top., 224 (2015), 1469–1476. https://doi.org/10.1140/epjst/e2015-02472-1 doi: 10.1140/epjst/e2015-02472-1
[21]	X. Wang, G. Chen, A chaotic system with only one stable equilibrium, Commun. Nonlinear Sci., 17 (2012), 1264–1272. https://doi.org/10.1016/j.cnsns.2011.07.017 doi: 10.1016/j.cnsns.2011.07.017
[22]	X. Wang, A. Akgul, S. Cicek, V. T. Pham, D. V. Hoang, A chaotic system with two stable equilibrium points: dynamics, circuit realization and communication application, Int. J. Bifurcat. Chaos, 27 (2017), 1750130. https://doi.org/10.1142/S0218127417501309 doi: 10.1142/S0218127417501309
[23]	P. C. Rech, Self-excited and hidden attractors in a multistable jerk system, Chaos Soliton. Fract., 164 (2022), 112614. https://doi.org/10.1016/j.chaos.2022.112614 doi: 10.1016/j.chaos.2022.112614
[24]	V. T. Pham, S. Jafari, C. Volos, T. Kapitaniak, A gallery of chaotic systems with an infinite number of equilibrium points, Chaos Soliton. Fract., 93 (2016), 58–63. https://doi.org/10.1016/j.chaos.2016.10.002 doi: 10.1016/j.chaos.2016.10.002
[25]	M. Molaie, S. Jafari, J. C. Sprott, S. M. Golpayegani, Simple chaotic flows with one stable equilibrium, Int. J. Bifurcat. Chaos, 23 (2013), 1350188. https://doi.org/10.1142/S0218127413501885 doi: 10.1142/S0218127413501885
[26]	S. Jafari, J. C. Sprott, Simple chaotic flows with a line equilibrium, Chaos Soliton. Fract., 57 (2013), 79–84. https://doi.org/10.1016/j.chaos.2013.08.018 doi: 10.1016/j.chaos.2013.08.018
[27]	S. Kumarasamy, M. Banerjee, V. Varshney, M. D. Shrimali, N. V. Kuznetsov, A. Prasad, Saddle-node bifurcation of periodic orbit route to hidden attractors, Phys. Rev. E, 107 (2023), L052201. https://doi.org/10.1103/PhysRevE.107.L052201 doi: 10.1103/PhysRevE.107.L052201
[28]	R. Balamurali, J. Kengne, R. G. Chengui, K. Rajagopal, Coupled van der Pol and Duffing oscillators: emergence of antimonotonicity and coexisting multiple self-excited and hidden oscillations, Eur. Phys. J. Plus, 137 (2022), 789. https://doi.org/10.1140/epjp/s13360-022-03000-2 doi: 10.1140/epjp/s13360-022-03000-2
[29]	B. Li, B. Sang, M. Liu, X. Hu, X. Zhang, N. Wang, Some jerk systems with hidden chaotic dynamics, Int. J. Bifurcat. Chaos, 33 (2023), 2350069. https://doi.org/10.1142/S0218127423500694 doi: 10.1142/S0218127423500694
[30]	C. Dong, M. Yang, L. Jia, Z. Li, Dynamics investigation and chaos-based application of a novel no-equilibrium system with coexisting hidden attractors, Phys. A, 633 (2024), 129391. https://doi.org/10.1016/j.physa.2023.129391 doi: 10.1016/j.physa.2023.129391
[31]	J. C. Sprott, Strange attractors with various equilibrium types, Eur. Phys. J. Spec. Top., 224 (2015), 1409–1419. https://doi.org/10.1140/epjst/e2015-02469-8 doi: 10.1140/epjst/e2015-02469-8
[32]	M. A. Jafari, E. Mliki, A. Akgul, V. T. Pham, S. T. Kingni, X. Wang, S. Jafari, Chameleon: the most hidden chaotic flow, Nonlinear Dyn., 88 (2017), 2303–2317. https://doi.org/10.1007/s11071-017-3378-4 doi: 10.1007/s11071-017-3378-4
[33]	K. Rajagopal, A. Karthikeyan, P. Duraisamy, Hyperchaotic chameleon: fractional order FPGA implementation, Complexity, 2017 (2017), 8979408. https://doi.org/10.1155/2017/8979408 doi: 10.1155/2017/8979408
[34]	H. Natiq, M. R. Said, M. R. Ariffin, S. He, L. Rondoni, S. Banerjee, Self-excited and hidden attractors in a novel chaotic system with complicated multistability, Eur. Phys. J. Plus, 133 (2018), 557. https://doi.org/10.1140/epjp/i2018-12360-y doi: 10.1140/epjp/i2018-12360-y
[35]	S. Cang, Y. Li, R. Zhang, Z. Wang, Hidden and self-excited coexisting attractors in a Lorenz-like system with two equilibrium points, Nonlinear Dyn., 95 (2019), 381–390. https://doi.org/10.1007/s11071-018-4570-x doi: 10.1007/s11071-018-4570-x
[36]	Q. Yang, Z. Wei, G. Chen, An unusual 3D autonomous quadratic chaotic system with two stable node-foci, Int. J. Bifurcat. Chaos, 20 (2010), 1061–1083. https://doi.org/10.1142/S0218127410026320 doi: 10.1142/S0218127410026320
[37]	V. F. Signing, G. G. Tegue, M. Kountchou, Z. T. Njitacke, N. Tsafack, J. D. Nkapkop, et al., A cryptosystem based on a chameleon chaotic system and dynamic DNA coding, Chaos Soliton. Fract., 155 (2022), 111777. https://doi.org/10.1016/j.chaos.2021.111777 doi: 10.1016/j.chaos.2021.111777
[38]	W. Fan, D. Xu, Z. Chen, N. Wang, Q. Xu, On two-parameter bifurcation and analog circuit implementation of a chameleon chaotic system, Phys. Scr., 99 (2023), 015218. https://doi.org/10.1088/1402-4896/ad1231 doi: 10.1088/1402-4896/ad1231
[39]	A. Tiwari, R. Nathasarma, B. K. Roy, A new time-reversible 3D chaotic system with coexisting dissipative and conservative behaviors and its active non-linear control, J. Franklin I., 361 (2024), 106637. https://doi.org/10.1016/j.jfranklin.2024.01.038 doi: 10.1016/j.jfranklin.2024.01.038
[40]	R. Zhou, Y. Gu, J. Cui, G. Ren, S. Yu, Nonlinear dynamic analysis of supercritical and subcritical Hopf bifurcations in gas foil bearing-rotor systems, Nonlinear Dyn., 103 (2021), 2241–2256. https://doi.org/10.1007/s11071-021-06234-4 doi: 10.1007/s11071-021-06234-4
[41]	N. V. Stankevich, N. V. Kuznetsov, G. A. Leonov, L. O. Chua, Scenario of the birth of hidden attractors in the Chua circuit, Int. J. Bifurcat. Chaos, 27 (2017), 1730038. https://doi.org/10.1142/S0218127417300385 doi: 10.1142/S0218127417300385
[42]	H. Zhao, Y. Lin, Y. Dai, Hopf bifurcation and hidden attractor of a modified Chua's equation, Nonlinear Dyn., 90 (2017), 2013–2021. https://doi.org/10.1007/s11071-017-3777-6 doi: 10.1007/s11071-017-3777-6
[43]	M. Liu, B. Sang, N. Wang, I. Ahmad, Chaotic dynamics by some quadratic jerk systems, Axioms, 10 (2021), 227. https://doi.org/10.3390/axioms10030227 doi: 10.3390/axioms10030227
[44]	Q. Yang, D. Zhu, L. Yang, A new 7D hyperchaotic system with five positive Lyapunov exponents coined, Int. J. Bifurcat. Chaos, 28 (2018), 1850057. https://doi.org/10.1142/S0218127418500578 doi: 10.1142/S0218127418500578
[45]	Z. Li, K. Chen, Neuromorphic behaviors in a neuron circuit based on current-controlled Chua Corsage Memristor, Chaos Soliton. Fract., 175 (2023), 114017. https://doi.org/10.1016/j.chaos.2023.114017 doi: 10.1016/j.chaos.2023.114017
[46]	Y. Liu, Y. Zhou, B. Guo, Hopf bifurcation, periodic solutions, and control of a new 4D hyperchaotic system, Mathematics, 11 (2023), 2699. https://doi.org/10.3390/math11122699 doi: 10.3390/math11122699
[47]	C. Li, J. C. Sprott, Variable-boostable chaotic flows, Optik, 127 (2016), 10389–10398. https://doi.org/10.1016/j.ijleo.2016.08.046 doi: 10.1016/j.ijleo.2016.08.046
[48]	C. Li, A. Akgul, L. Bi, Y. Xu, C. Zhang, A chaotic jerk oscillator with interlocked offset boosting, Eur. Phys. J. Plus, 139 (2024), 242. https://doi.org/10.1140/epjp/s13360-024-05040-2 doi: 10.1140/epjp/s13360-024-05040-2
[49]	C. Li, X. Wang, G. Chen, Diagnosing multistability by offset boosting, Nonlinear Dyn., 90 (2017), 1335–1341. https://doi.org/10.1007/s11071-017-3729-1 doi: 10.1007/s11071-017-3729-1
[50]	X. Gao, Hamilton energy of a complex chaotic system and offset boosting, Phys. Scr., 99 (2024), 015244. https://doi.org/10.1088/1402-4896/ad1739 doi: 10.1088/1402-4896/ad1739
[51]	X. Zhang, C. Li, T. Lei, H. Fu, Z. Liu, Offset boosting in a memristive hyperchaotic system, Phys. Scr., 99 (2024), 015247. https://doi.org/10.1088/1402-4896/ad156e doi: 10.1088/1402-4896/ad156e
[52]	H. Dang-Vu, C. Delcarte, Bifurcations et Chaos: une introduction à la dynamique contemporaine avec des programmes en Pascal, Fortran et Mathematica, Paris: Ellipses, 2000.
[53]	C. Grebogi, E. Ott, J. A. Yorke, Chaos, strange attractors, and fractal basin boundaries in nonlinear dynamics, Science, 238 (1987), 632–638. https://doi.org/10.1126/science.238.4827.632 doi: 10.1126/science.238.4827.632

This article has been cited by:

1.	Amin Ghoreyshi, Mostafa Abbaszadeh, Mahmoud A. Zaky, Mehdi Dehghan, Finite block method for nonlinear time-fractional partial integro-differential equations: Stability, convergence, and numerical analysis, 2025, 214, 01689274, 82, 10.1016/j.apnum.2025.03.002
2.	Mostafa Abbaszadeh, Mahmoud A. Zaky, Mehdi Dehghan, Virtual element approximation and BDF2 time-discrete scheme for a partial integro-differential equation with a singular Abel's kernel, 2025, 501, 00963003, 129451, 10.1016/j.amc.2025.129451
3.	Ali H. Tedjani, Sharifah E. Alhazmi, Samer S. Ezz-Eldien, An operational approach for one- and two-dimension high-order multi-pantograph Volterra integro-differential equation, 2025, 10, 2473-6988, 9274, 10.3934/math.2025426

Reader Comments

Your name:*

Email:*
© 2025 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)